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 whose 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:   PetscFunctionBegin;
 31:   while (len--) {
 32:     max -= tup[len];
 33:     if (!max) {
 34:       tup[len] = 0;
 35:       break;
 36:     }
 37:   }
 38:   tup[++len]++;
 39:   PetscFunctionReturn(PETSC_SUCCESS);
 40: }

 42: /*
 43:   PetscDualSpaceTensorPointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that are all less than or equal to 'max'.
 44:                                                     Ordering is lexicographic with lowest index as least significant in ordering.
 45:                                                     e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {2,1}, {0,2}, {1,2}, {2,2}.

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

 52:   Output Parameter:
 53: . tup - A tuple of `len` integers whose entries are at most `max`

 55:   Level: developer

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

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

 75: /*@C
 76:   PetscDualSpaceRegister - Adds a new `PetscDualSpaceType`

 78:   Not Collective

 80:   Input Parameters:
 81: + sname    - The name of a new user-defined creation routine
 82: - function - The creation routine

 84:   Example Usage:
 85: .vb
 86:     PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
 87: .ve

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

 99:   Level: advanced

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

104: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceRegisterAll()`, `PetscDualSpaceRegisterDestroy()`
105: @*/
106: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
107: {
108:   PetscFunctionBegin;
109:   PetscCall(PetscFunctionListAdd(&PetscDualSpaceList, sname, function));
110:   PetscFunctionReturn(PETSC_SUCCESS);
111: }

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

116:   Collective

118:   Input Parameters:
119: + sp   - The `PetscDualSpace` object
120: - name - The kind of space

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

125:   Level: intermediate

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

134:   PetscFunctionBegin;
136:   PetscCall(PetscObjectTypeCompare((PetscObject)sp, name, &match));
137:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

139:   if (!PetscDualSpaceRegisterAllCalled) PetscCall(PetscDualSpaceRegisterAll());
140:   PetscCall(PetscFunctionListFind(PetscDualSpaceList, name, &r));
141:   PetscCheck(r, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDualSpace type: %s", name);

143:   PetscTryTypeMethod(sp, destroy);
144:   sp->ops->destroy = NULL;

146:   PetscCall((*r)(sp));
147:   PetscCall(PetscObjectChangeTypeName((PetscObject)sp, name));
148:   PetscFunctionReturn(PETSC_SUCCESS);
149: }

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

154:   Not Collective

156:   Input Parameter:
157: . sp - The `PetscDualSpace`

159:   Output Parameter:
160: . name - The `PetscDualSpaceType` name

162:   Level: intermediate

164: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceSetType()`, `PetscDualSpaceCreate()`
165: @*/
166: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
167: {
168:   PetscFunctionBegin;
170:   PetscAssertPointer(name, 2);
171:   if (!PetscDualSpaceRegisterAllCalled) PetscCall(PetscDualSpaceRegisterAll());
172:   *name = ((PetscObject)sp)->type_name;
173:   PetscFunctionReturn(PETSC_SUCCESS);
174: }

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

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

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

209:   Collective

211:   Input Parameters:
212: + A    - the `PetscDualSpace` object
213: . obj  - Optional object, provides the options prefix
214: - name - command line option name

216:   Level: intermediate

218:   Note:
219:   See `PetscObjectViewFromOptions()` for possible command line values

221: .seealso: `PetscDualSpace`, `PetscDualSpaceView()`, `PetscObjectViewFromOptions()`, `PetscDualSpaceCreate()`
222: @*/
223: PetscErrorCode PetscDualSpaceViewFromOptions(PetscDualSpace A, PetscObject obj, const char name[])
224: {
225:   PetscFunctionBegin;
227:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
228:   PetscFunctionReturn(PETSC_SUCCESS);
229: }

231: /*@
232:   PetscDualSpaceView - Views a `PetscDualSpace`

234:   Collective

236:   Input Parameters:
237: + sp - the `PetscDualSpace` object to view
238: - v  - the viewer

240:   Level: beginner

242: .seealso: `PetscViewer`, `PetscDualSpaceDestroy()`, `PetscDualSpace`
243: @*/
244: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
245: {
246:   PetscBool iascii;

248:   PetscFunctionBegin;
251:   if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &v));
252:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii));
253:   if (iascii) PetscCall(PetscDualSpaceView_ASCII(sp, v));
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

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

260:   Collective

262:   Input Parameter:
263: . sp - the `PetscDualSpace` object to set options for

265:   Options Database Keys:
266: + -petscdualspace_order <order>                 - the approximation order of the space
267: . -petscdualspace_form_degree <deg>             - the form degree, say 0 for point evaluations, or 2 for area integrals
268: . -petscdualspace_components <c>                - the number of components, say d for a vector field
269: . -petscdualspace_refcell <celltype>            - Reference cell type name
270: . -petscdualspace_lagrange_continuity           - Flag for continuous element
271: . -petscdualspace_lagrange_tensor               - Flag for tensor dual space
272: . -petscdualspace_lagrange_trimmed              - Flag for trimmed dual space
273: . -petscdualspace_lagrange_node_type <nodetype> - Lagrange node location type
274: . -petscdualspace_lagrange_node_endpoints       - Flag for nodes that include endpoints
275: . -petscdualspace_lagrange_node_exponent        - Gauss-Jacobi weight function exponent
276: . -petscdualspace_lagrange_use_moments          - Use moments (where appropriate) for functionals
277: - -petscdualspace_lagrange_moment_order <order> - Quadrature order for moment functionals

279:   Level: intermediate

281: .seealso: `PetscDualSpaceView()`, `PetscDualSpace`, `PetscObjectSetFromOptions()`
282: @*/
283: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
284: {
285:   DMPolytopeType refCell = DM_POLYTOPE_TRIANGLE;
286:   const char    *defaultType;
287:   char           name[256];
288:   PetscBool      flg;

290:   PetscFunctionBegin;
292:   if (!((PetscObject)sp)->type_name) {
293:     defaultType = PETSCDUALSPACELAGRANGE;
294:   } else {
295:     defaultType = ((PetscObject)sp)->type_name;
296:   }
297:   if (!PetscSpaceRegisterAllCalled) PetscCall(PetscSpaceRegisterAll());

299:   PetscObjectOptionsBegin((PetscObject)sp);
300:   PetscCall(PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg));
301:   if (flg) {
302:     PetscCall(PetscDualSpaceSetType(sp, name));
303:   } else if (!((PetscObject)sp)->type_name) {
304:     PetscCall(PetscDualSpaceSetType(sp, defaultType));
305:   }
306:   PetscCall(PetscOptionsBoundedInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL, 0));
307:   PetscCall(PetscOptionsInt("-petscdualspace_form_degree", "The form degree of the dofs", "PetscDualSpaceSetFormDegree", sp->k, &sp->k, NULL));
308:   PetscCall(PetscOptionsBoundedInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL, 1));
309:   PetscTryTypeMethod(sp, setfromoptions, PetscOptionsObject);
310:   PetscCall(PetscOptionsEnum("-petscdualspace_refcell", "Reference cell shape", "PetscDualSpaceSetReferenceCell", DMPolytopeTypes, (PetscEnum)refCell, (PetscEnum *)&refCell, &flg));
311:   if (flg) {
312:     DM K;

314:     PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, refCell, &K));
315:     PetscCall(PetscDualSpaceSetDM(sp, K));
316:     PetscCall(DMDestroy(&K));
317:   }

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

326: /*@
327:   PetscDualSpaceSetUp - Construct a basis for a `PetscDualSpace`

329:   Collective

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

334:   Level: intermediate

336: .seealso: `PetscDualSpaceView()`, `PetscDualSpaceDestroy()`, `PetscDualSpace`
337: @*/
338: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
339: {
340:   PetscFunctionBegin;
342:   if (sp->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
343:   PetscCall(PetscLogEventBegin(PETSCDUALSPACE_SetUp, sp, 0, 0, 0));
344:   sp->setupcalled = PETSC_TRUE;
345:   PetscTryTypeMethod(sp, setup);
346:   PetscCall(PetscLogEventEnd(PETSCDUALSPACE_SetUp, sp, 0, 0, 0));
347:   if (sp->setfromoptionscalled) PetscCall(PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view"));
348:   PetscFunctionReturn(PETSC_SUCCESS);
349: }

351: static PetscErrorCode PetscDualSpaceClearDMData_Internal(PetscDualSpace sp, DM dm)
352: {
353:   PetscInt pStart = -1, pEnd = -1, depth = -1;

355:   PetscFunctionBegin;
356:   if (!dm) PetscFunctionReturn(PETSC_SUCCESS);
357:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
358:   PetscCall(DMPlexGetDepth(dm, &depth));

360:   if (sp->pointSpaces) {
361:     PetscInt i;

363:     for (i = 0; i < pEnd - pStart; i++) PetscCall(PetscDualSpaceDestroy(&(sp->pointSpaces[i])));
364:   }
365:   PetscCall(PetscFree(sp->pointSpaces));

367:   if (sp->heightSpaces) {
368:     PetscInt i;

370:     for (i = 0; i <= depth; i++) PetscCall(PetscDualSpaceDestroy(&(sp->heightSpaces[i])));
371:   }
372:   PetscCall(PetscFree(sp->heightSpaces));

374:   PetscCall(PetscSectionDestroy(&(sp->pointSection)));
375:   PetscCall(PetscSectionDestroy(&(sp->intPointSection)));
376:   PetscCall(PetscQuadratureDestroy(&(sp->intNodes)));
377:   PetscCall(VecDestroy(&(sp->intDofValues)));
378:   PetscCall(VecDestroy(&(sp->intNodeValues)));
379:   PetscCall(MatDestroy(&(sp->intMat)));
380:   PetscCall(PetscQuadratureDestroy(&(sp->allNodes)));
381:   PetscCall(VecDestroy(&(sp->allDofValues)));
382:   PetscCall(VecDestroy(&(sp->allNodeValues)));
383:   PetscCall(MatDestroy(&(sp->allMat)));
384:   PetscCall(PetscFree(sp->numDof));
385:   PetscFunctionReturn(PETSC_SUCCESS);
386: }

388: /*@
389:   PetscDualSpaceDestroy - Destroys a `PetscDualSpace` object

391:   Collective

393:   Input Parameter:
394: . sp - the `PetscDualSpace` object to destroy

396:   Level: beginner

398: .seealso: `PetscDualSpace`, `PetscDualSpaceView()`, `PetscDualSpace()`, `PetscDualSpaceCreate()`
399: @*/
400: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
401: {
402:   PetscInt dim, f;
403:   DM       dm;

405:   PetscFunctionBegin;
406:   if (!*sp) PetscFunctionReturn(PETSC_SUCCESS);

409:   if (--((PetscObject)(*sp))->refct > 0) {
410:     *sp = NULL;
411:     PetscFunctionReturn(PETSC_SUCCESS);
412:   }
413:   ((PetscObject)(*sp))->refct = 0;

415:   PetscCall(PetscDualSpaceGetDimension(*sp, &dim));
416:   dm = (*sp)->dm;

418:   PetscTryTypeMethod((*sp), destroy);
419:   PetscCall(PetscDualSpaceClearDMData_Internal(*sp, dm));

421:   for (f = 0; f < dim; ++f) PetscCall(PetscQuadratureDestroy(&(*sp)->functional[f]));
422:   PetscCall(PetscFree((*sp)->functional));
423:   PetscCall(DMDestroy(&(*sp)->dm));
424:   PetscCall(PetscHeaderDestroy(sp));
425:   PetscFunctionReturn(PETSC_SUCCESS);
426: }

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

431:   Collective

433:   Input Parameter:
434: . comm - The communicator for the `PetscDualSpace` object

436:   Output Parameter:
437: . sp - The `PetscDualSpace` object

439:   Level: beginner

441: .seealso: `PetscDualSpace`, `PetscDualSpaceSetType()`, `PETSCDUALSPACELAGRANGE`
442: @*/
443: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
444: {
445:   PetscDualSpace s;

447:   PetscFunctionBegin;
448:   PetscAssertPointer(sp, 2);
449:   PetscCall(PetscCitationsRegister(FECitation, &FEcite));
450:   *sp = NULL;
451:   PetscCall(PetscFEInitializePackage());

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

455:   s->order       = 0;
456:   s->Nc          = 1;
457:   s->k           = 0;
458:   s->spdim       = -1;
459:   s->spintdim    = -1;
460:   s->uniform     = PETSC_TRUE;
461:   s->setupcalled = PETSC_FALSE;

463:   *sp = s;
464:   PetscFunctionReturn(PETSC_SUCCESS);
465: }

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

470:   Collective

472:   Input Parameter:
473: . sp - The original `PetscDualSpace`

475:   Output Parameter:
476: . spNew - The duplicate `PetscDualSpace`

478:   Level: beginner

480: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
481: @*/
482: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
483: {
484:   DM                 dm;
485:   PetscDualSpaceType type;
486:   const char        *name;

488:   PetscFunctionBegin;
490:   PetscAssertPointer(spNew, 2);
491:   PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)sp), spNew));
492:   name = ((PetscObject)sp)->name;
493:   if (name) { PetscCall(PetscObjectSetName((PetscObject)*spNew, name)); }
494:   PetscCall(PetscDualSpaceGetType(sp, &type));
495:   PetscCall(PetscDualSpaceSetType(*spNew, type));
496:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
497:   PetscCall(PetscDualSpaceSetDM(*spNew, dm));

499:   (*spNew)->order   = sp->order;
500:   (*spNew)->k       = sp->k;
501:   (*spNew)->Nc      = sp->Nc;
502:   (*spNew)->uniform = sp->uniform;
503:   PetscTryTypeMethod(sp, duplicate, *spNew);
504:   PetscFunctionReturn(PETSC_SUCCESS);
505: }

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

510:   Not Collective

512:   Input Parameter:
513: . sp - The `PetscDualSpace`

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

518:   Level: intermediate

520: .seealso: `PetscDualSpace`, `PetscDualSpaceSetDM()`, `PetscDualSpaceCreate()`
521: @*/
522: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
523: {
524:   PetscFunctionBegin;
526:   PetscAssertPointer(dm, 2);
527:   *dm = sp->dm;
528:   PetscFunctionReturn(PETSC_SUCCESS);
529: }

531: /*@
532:   PetscDualSpaceSetDM - Get the `DM` representing the reference cell

534:   Not Collective

536:   Input Parameters:
537: + sp - The `PetscDual`Space
538: - dm - The reference cell

540:   Level: intermediate

542: .seealso: `PetscDualSpace`, `DM`, `PetscDualSpaceGetDM()`, `PetscDualSpaceCreate()`
543: @*/
544: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
545: {
546:   PetscFunctionBegin;
549:   PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change DM after dualspace is set up");
550:   PetscCall(PetscObjectReference((PetscObject)dm));
551:   if (sp->dm && sp->dm != dm) PetscCall(PetscDualSpaceClearDMData_Internal(sp, sp->dm));
552:   PetscCall(DMDestroy(&sp->dm));
553:   sp->dm = dm;
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }

557: /*@
558:   PetscDualSpaceGetOrder - Get the order of the dual space

560:   Not Collective

562:   Input Parameter:
563: . sp - The `PetscDualSpace`

565:   Output Parameter:
566: . order - The order

568:   Level: intermediate

570: .seealso: `PetscDualSpace`, `PetscDualSpaceSetOrder()`, `PetscDualSpaceCreate()`
571: @*/
572: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
573: {
574:   PetscFunctionBegin;
576:   PetscAssertPointer(order, 2);
577:   *order = sp->order;
578:   PetscFunctionReturn(PETSC_SUCCESS);
579: }

581: /*@
582:   PetscDualSpaceSetOrder - Set the order of the dual space

584:   Not Collective

586:   Input Parameters:
587: + sp    - The `PetscDualSpace`
588: - order - The order

590:   Level: intermediate

592: .seealso: `PetscDualSpace`, `PetscDualSpaceGetOrder()`, `PetscDualSpaceCreate()`
593: @*/
594: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
595: {
596:   PetscFunctionBegin;
598:   PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change order after dualspace is set up");
599:   sp->order = order;
600:   PetscFunctionReturn(PETSC_SUCCESS);
601: }

603: /*@
604:   PetscDualSpaceGetNumComponents - Return the number of components for this space

606:   Input Parameter:
607: . sp - The `PetscDualSpace`

609:   Output Parameter:
610: . Nc - The number of components

612:   Level: intermediate

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

617: .seealso: `PetscDualSpaceSetNumComponents()`, `PetscDualSpaceGetDimension()`, `PetscDualSpaceCreate()`, `PetscDualSpace`
618: @*/
619: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
620: {
621:   PetscFunctionBegin;
623:   PetscAssertPointer(Nc, 2);
624:   *Nc = sp->Nc;
625:   PetscFunctionReturn(PETSC_SUCCESS);
626: }

628: /*@
629:   PetscDualSpaceSetNumComponents - Set the number of components for this space

631:   Input Parameters:
632: + sp - The `PetscDualSpace`
633: - Nc - The number of components

635:   Level: intermediate

637: .seealso: `PetscDualSpaceGetNumComponents()`, `PetscDualSpaceCreate()`, `PetscDualSpace`
638: @*/
639: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
640: {
641:   PetscFunctionBegin;
643:   PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up");
644:   sp->Nc = Nc;
645:   PetscFunctionReturn(PETSC_SUCCESS);
646: }

648: /*@
649:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

651:   Not Collective

653:   Input Parameters:
654: + sp - The `PetscDualSpace`
655: - i  - The basis number

657:   Output Parameter:
658: . functional - The basis functional

660:   Level: intermediate

662: .seealso: `PetscDualSpace`, `PetscQuadrature`, `PetscDualSpaceGetDimension()`, `PetscDualSpaceCreate()`
663: @*/
664: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
665: {
666:   PetscInt dim;

668:   PetscFunctionBegin;
670:   PetscAssertPointer(functional, 3);
671:   PetscCall(PetscDualSpaceGetDimension(sp, &dim));
672:   PetscCheck(!(i < 0) && !(i >= dim), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", i, dim);
673:   *functional = sp->functional[i];
674:   PetscFunctionReturn(PETSC_SUCCESS);
675: }

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

680:   Not Collective

682:   Input Parameter:
683: . sp - The `PetscDualSpace`

685:   Output Parameter:
686: . dim - The dimension

688:   Level: intermediate

690: .seealso: `PetscDualSpace`, `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
691: @*/
692: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
693: {
694:   PetscFunctionBegin;
696:   PetscAssertPointer(dim, 2);
697:   if (sp->spdim < 0) {
698:     PetscSection section;

700:     PetscCall(PetscDualSpaceGetSection(sp, &section));
701:     if (section) {
702:       PetscCall(PetscSectionGetStorageSize(section, &(sp->spdim)));
703:     } else sp->spdim = 0;
704:   }
705:   *dim = sp->spdim;
706:   PetscFunctionReturn(PETSC_SUCCESS);
707: }

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

712:   Not Collective

714:   Input Parameter:
715: . sp - The `PetscDualSpace`

717:   Output Parameter:
718: . intdim - The dimension

720:   Level: intermediate

722: .seealso: `PetscDualSpace`, `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
723: @*/
724: PetscErrorCode PetscDualSpaceGetInteriorDimension(PetscDualSpace sp, PetscInt *intdim)
725: {
726:   PetscFunctionBegin;
728:   PetscAssertPointer(intdim, 2);
729:   if (sp->spintdim < 0) {
730:     PetscSection section;

732:     PetscCall(PetscDualSpaceGetSection(sp, &section));
733:     if (section) {
734:       PetscCall(PetscSectionGetConstrainedStorageSize(section, &(sp->spintdim)));
735:     } else sp->spintdim = 0;
736:   }
737:   *intdim = sp->spintdim;
738:   PetscFunctionReturn(PETSC_SUCCESS);
739: }

741: /*@
742:   PetscDualSpaceGetUniform - Whether this dual space is uniform

744:   Not Collective

746:   Input Parameter:
747: . sp - A dual space

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

753:   Level: advanced

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

759: .seealso: `PetscDualSpace`, `PetscDualSpaceGetPointSubspace()`, `PetscDualSpaceGetSymmetries()`
760: @*/
761: PetscErrorCode PetscDualSpaceGetUniform(PetscDualSpace sp, PetscBool *uniform)
762: {
763:   PetscFunctionBegin;
765:   PetscAssertPointer(uniform, 2);
766:   *uniform = sp->uniform;
767:   PetscFunctionReturn(PETSC_SUCCESS);
768: }

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

773:   Not Collective

775:   Input Parameter:
776: . sp - The `PetscDualSpace`

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

781:   Level: intermediate

783: .seealso: `PetscDualSpace`, `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
784: @*/
785: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
786: {
787:   PetscFunctionBegin;
789:   PetscAssertPointer(numDof, 2);
790:   PetscCheck(sp->uniform, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "A non-uniform space does not have a fixed number of dofs for each height");
791:   if (!sp->numDof) {
792:     DM           dm;
793:     PetscInt     depth, d;
794:     PetscSection section;

796:     PetscCall(PetscDualSpaceGetDM(sp, &dm));
797:     PetscCall(DMPlexGetDepth(dm, &depth));
798:     PetscCall(PetscCalloc1(depth + 1, &(sp->numDof)));
799:     PetscCall(PetscDualSpaceGetSection(sp, &section));
800:     for (d = 0; d <= depth; d++) {
801:       PetscInt dStart, dEnd;

803:       PetscCall(DMPlexGetDepthStratum(dm, d, &dStart, &dEnd));
804:       if (dEnd <= dStart) continue;
805:       PetscCall(PetscSectionGetDof(section, dStart, &(sp->numDof[d])));
806:     }
807:   }
808:   *numDof = sp->numDof;
809:   PetscCheck(*numDof, PetscObjectComm((PetscObject)sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
810:   PetscFunctionReturn(PETSC_SUCCESS);
811: }

813: /* create the section of the right size and set a permutation for topological ordering */
814: PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace sp, PetscSection *topSection)
815: {
816:   DM           dm;
817:   PetscInt     pStart, pEnd, cStart, cEnd, c, depth, count, i;
818:   PetscInt    *seen, *perm;
819:   PetscSection section;

821:   PetscFunctionBegin;
822:   dm = sp->dm;
823:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &section));
824:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
825:   PetscCall(PetscSectionSetChart(section, pStart, pEnd));
826:   PetscCall(PetscCalloc1(pEnd - pStart, &seen));
827:   PetscCall(PetscMalloc1(pEnd - pStart, &perm));
828:   PetscCall(DMPlexGetDepth(dm, &depth));
829:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
830:   for (c = cStart, count = 0; c < cEnd; c++) {
831:     PetscInt  closureSize = -1, e;
832:     PetscInt *closure     = NULL;

834:     perm[count++]    = c;
835:     seen[c - pStart] = 1;
836:     PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
837:     for (e = 0; e < closureSize; e++) {
838:       PetscInt point = closure[2 * e];

840:       if (seen[point - pStart]) continue;
841:       perm[count++]        = point;
842:       seen[point - pStart] = 1;
843:     }
844:     PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
845:   }
846:   PetscCheck(count == pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad topological ordering");
847:   for (i = 0; i < pEnd - pStart; i++)
848:     if (perm[i] != i) break;
849:   if (i < pEnd - pStart) {
850:     IS permIS;

852:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
853:     PetscCall(ISSetPermutation(permIS));
854:     PetscCall(PetscSectionSetPermutation(section, permIS));
855:     PetscCall(ISDestroy(&permIS));
856:   } else {
857:     PetscCall(PetscFree(perm));
858:   }
859:   PetscCall(PetscFree(seen));
860:   *topSection = section;
861:   PetscFunctionReturn(PETSC_SUCCESS);
862: }

864: /* mark boundary points and set up */
865: PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace sp, PetscSection section)
866: {
867:   DM       dm;
868:   DMLabel  boundary;
869:   PetscInt pStart, pEnd, p;

871:   PetscFunctionBegin;
872:   dm = sp->dm;
873:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "boundary", &boundary));
874:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
875:   PetscCall(DMPlexMarkBoundaryFaces(dm, 1, boundary));
876:   PetscCall(DMPlexLabelComplete(dm, boundary));
877:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
878:   for (p = pStart; p < pEnd; p++) {
879:     PetscInt bval;

881:     PetscCall(DMLabelGetValue(boundary, p, &bval));
882:     if (bval == 1) {
883:       PetscInt dof;

885:       PetscCall(PetscSectionGetDof(section, p, &dof));
886:       PetscCall(PetscSectionSetConstraintDof(section, p, dof));
887:     }
888:   }
889:   PetscCall(DMLabelDestroy(&boundary));
890:   PetscCall(PetscSectionSetUp(section));
891:   PetscFunctionReturn(PETSC_SUCCESS);
892: }

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

897:   Collective

899:   Input Parameter:
900: . sp - The `PetscDualSpace`

902:   Output Parameter:
903: . section - The section

905:   Level: advanced

907: .seealso: `PetscDualSpace`, `PetscSection`, `PetscDualSpaceCreate()`, `DMPLEX`
908: @*/
909: PetscErrorCode PetscDualSpaceGetSection(PetscDualSpace sp, PetscSection *section)
910: {
911:   PetscInt pStart, pEnd, p;

913:   PetscFunctionBegin;
914:   if (!sp->dm) {
915:     *section = NULL;
916:     PetscFunctionReturn(PETSC_SUCCESS);
917:   }
918:   if (!sp->pointSection) {
919:     /* mark the boundary */
920:     PetscCall(PetscDualSpaceSectionCreate_Internal(sp, &(sp->pointSection)));
921:     PetscCall(DMPlexGetChart(sp->dm, &pStart, &pEnd));
922:     for (p = pStart; p < pEnd; p++) {
923:       PetscDualSpace psp;

925:       PetscCall(PetscDualSpaceGetPointSubspace(sp, p, &psp));
926:       if (psp) {
927:         PetscInt dof;

929:         PetscCall(PetscDualSpaceGetInteriorDimension(psp, &dof));
930:         PetscCall(PetscSectionSetDof(sp->pointSection, p, dof));
931:       }
932:     }
933:     PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, sp->pointSection));
934:   }
935:   *section = sp->pointSection;
936:   PetscFunctionReturn(PETSC_SUCCESS);
937: }

939: /*@
940:   PetscDualSpaceGetInteriorSection - Create a `PetscSection` over the reference cell with the layout from this space
941:   for interior degrees of freedom

943:   Collective

945:   Input Parameter:
946: . sp - The `PetscDualSpace`

948:   Output Parameter:
949: . section - The interior section

951:   Level: advanced

953:   Note:
954:   Most reference domains have one cell, in which case the only cell will have
955:   all of the interior degrees of freedom in the interior section.  But
956:   for `PETSCDUALSPACEREFINED` there may be other mesh points in the interior,
957:   and this section describes their layout.

959: .seealso: `PetscDualSpace`, `PetscSection`, `PetscDualSpaceCreate()`, `DMPLEX`
960: @*/
961: PetscErrorCode PetscDualSpaceGetInteriorSection(PetscDualSpace sp, PetscSection *section)
962: {
963:   PetscInt pStart, pEnd, p;

965:   PetscFunctionBegin;
966:   if (!sp->dm) {
967:     *section = NULL;
968:     PetscFunctionReturn(PETSC_SUCCESS);
969:   }
970:   if (!sp->intPointSection) {
971:     PetscSection full_section;

973:     PetscCall(PetscDualSpaceGetSection(sp, &full_section));
974:     PetscCall(PetscDualSpaceSectionCreate_Internal(sp, &(sp->intPointSection)));
975:     PetscCall(PetscSectionGetChart(full_section, &pStart, &pEnd));
976:     for (p = pStart; p < pEnd; p++) {
977:       PetscInt dof, cdof;

979:       PetscCall(PetscSectionGetDof(full_section, p, &dof));
980:       PetscCall(PetscSectionGetConstraintDof(full_section, p, &cdof));
981:       PetscCall(PetscSectionSetDof(sp->intPointSection, p, dof - cdof));
982:     }
983:     PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, sp->intPointSection));
984:   }
985:   *section = sp->intPointSection;
986:   PetscFunctionReturn(PETSC_SUCCESS);
987: }

989: /* this assumes that all of the point dual spaces store their interior dofs first, which is true when the point DMs
990:  * have one cell */
991: PetscErrorCode PetscDualSpacePushForwardSubspaces_Internal(PetscDualSpace sp, PetscInt sStart, PetscInt sEnd)
992: {
993:   PetscReal   *sv0, *v0, *J;
994:   PetscSection section;
995:   PetscInt     dim, s, k;
996:   DM           dm;

998:   PetscFunctionBegin;
999:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
1000:   PetscCall(DMGetDimension(dm, &dim));
1001:   PetscCall(PetscDualSpaceGetSection(sp, &section));
1002:   PetscCall(PetscMalloc3(dim, &v0, dim, &sv0, dim * dim, &J));
1003:   PetscCall(PetscDualSpaceGetFormDegree(sp, &k));
1004:   for (s = sStart; s < sEnd; s++) {
1005:     PetscReal      detJ, hdetJ;
1006:     PetscDualSpace ssp;
1007:     PetscInt       dof, off, f, sdim;
1008:     PetscInt       i, j;
1009:     DM             sdm;

1011:     PetscCall(PetscDualSpaceGetPointSubspace(sp, s, &ssp));
1012:     if (!ssp) continue;
1013:     PetscCall(PetscSectionGetDof(section, s, &dof));
1014:     PetscCall(PetscSectionGetOffset(section, s, &off));
1015:     /* get the first vertex of the reference cell */
1016:     PetscCall(PetscDualSpaceGetDM(ssp, &sdm));
1017:     PetscCall(DMGetDimension(sdm, &sdim));
1018:     PetscCall(DMPlexComputeCellGeometryAffineFEM(sdm, 0, sv0, NULL, NULL, &hdetJ));
1019:     PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, s, v0, J, NULL, &detJ));
1020:     /* compactify Jacobian */
1021:     for (i = 0; i < dim; i++)
1022:       for (j = 0; j < sdim; j++) J[i * sdim + j] = J[i * dim + j];
1023:     for (f = 0; f < dof; f++) {
1024:       PetscQuadrature fn;

1026:       PetscCall(PetscDualSpaceGetFunctional(ssp, f, &fn));
1027:       PetscCall(PetscQuadraturePushForward(fn, dim, sv0, v0, J, k, &(sp->functional[off + f])));
1028:     }
1029:   }
1030:   PetscCall(PetscFree3(v0, sv0, J));
1031:   PetscFunctionReturn(PETSC_SUCCESS);
1032: }

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

1037:   Input Parameters:
1038: + sp      - The `PetscDualSpace` object
1039: . f       - The basis functional index
1040: . time    - The time
1041: . 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)
1042: . numComp - The number of components for the function
1043: . func    - The input function
1044: - ctx     - A context for the function

1046:   Output Parameter:
1047: . value - numComp output values

1049:   Calling sequence:
1050: .vb
1051:   PetscErrorCode func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt numComponents, PetscScalar values[], void *ctx)
1052: .ve

1054:   Level: beginner

1056: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1057: @*/
1058: 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)
1059: {
1060:   PetscFunctionBegin;
1062:   PetscAssertPointer(cgeom, 4);
1063:   PetscAssertPointer(value, 8);
1064:   PetscUseTypeMethod(sp, apply, f, time, cgeom, numComp, func, ctx, value);
1065:   PetscFunctionReturn(PETSC_SUCCESS);
1066: }

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

1071:   Input Parameters:
1072: + sp        - The `PetscDualSpace` object
1073: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetAllData()`

1075:   Output Parameter:
1076: . spValue - The values of all dual space functionals

1078:   Level: advanced

1080: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1081: @*/
1082: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1083: {
1084:   PetscFunctionBegin;
1086:   PetscUseTypeMethod(sp, applyall, pointEval, spValue);
1087:   PetscFunctionReturn(PETSC_SUCCESS);
1088: }

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

1093:   Input Parameters:
1094: + sp        - The `PetscDualSpace` object
1095: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetInteriorData()`

1097:   Output Parameter:
1098: . spValue - The values of interior dual space functionals

1100:   Level: advanced

1102: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1103: @*/
1104: PetscErrorCode PetscDualSpaceApplyInterior(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1105: {
1106:   PetscFunctionBegin;
1108:   PetscUseTypeMethod(sp, applyint, pointEval, spValue);
1109:   PetscFunctionReturn(PETSC_SUCCESS);
1110: }

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

1115:   Input Parameters:
1116: + sp    - The `PetscDualSpace` object
1117: . f     - The basis functional index
1118: . time  - The time
1119: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1120: . Nc    - The number of components for the function
1121: . func  - The input function
1122: - ctx   - A context for the function

1124:   Output Parameter:
1125: . value - The output value

1127:   Calling sequence:
1128: .vb
1129:    PetscErrorCode func(PetscInt dim, PetscReal time, const PetscReal x[],PetscInt numComponents, PetscScalar values[], void *ctx)
1130: .ve

1132:   Level: advanced

1134:   Note:
1135:   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.

1137: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1138: @*/
1139: 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)
1140: {
1141:   DM               dm;
1142:   PetscQuadrature  n;
1143:   const PetscReal *points, *weights;
1144:   PetscReal        x[3];
1145:   PetscScalar     *val;
1146:   PetscInt         dim, dE, qNc, c, Nq, q;
1147:   PetscBool        isAffine;

1149:   PetscFunctionBegin;
1151:   PetscAssertPointer(value, 8);
1152:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
1153:   PetscCall(PetscDualSpaceGetFunctional(sp, f, &n));
1154:   PetscCall(PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights));
1155:   PetscCheck(dim == cgeom->dim, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %" PetscInt_FMT " != cell geometry dimension %" PetscInt_FMT, dim, cgeom->dim);
1156:   PetscCheck(qNc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_SIZ, "The quadrature components %" PetscInt_FMT " != function components %" PetscInt_FMT, qNc, Nc);
1157:   PetscCall(DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val));
1158:   *value   = 0.0;
1159:   isAffine = cgeom->isAffine;
1160:   dE       = cgeom->dimEmbed;
1161:   for (q = 0; q < Nq; ++q) {
1162:     if (isAffine) {
1163:       CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q * dim], x);
1164:       PetscCall((*func)(dE, time, x, Nc, val, ctx));
1165:     } else {
1166:       PetscCall((*func)(dE, time, &cgeom->v[dE * q], Nc, val, ctx));
1167:     }
1168:     for (c = 0; c < Nc; ++c) *value += val[c] * weights[q * Nc + c];
1169:   }
1170:   PetscCall(DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val));
1171:   PetscFunctionReturn(PETSC_SUCCESS);
1172: }

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

1177:   Input Parameters:
1178: + sp        - The `PetscDualSpace` object
1179: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetAllData()`

1181:   Output Parameter:
1182: . spValue - The values of all dual space functionals

1184:   Level: advanced

1186: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1187: @*/
1188: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1189: {
1190:   Vec pointValues, dofValues;
1191:   Mat allMat;

1193:   PetscFunctionBegin;
1195:   PetscAssertPointer(pointEval, 2);
1196:   PetscAssertPointer(spValue, 3);
1197:   PetscCall(PetscDualSpaceGetAllData(sp, NULL, &allMat));
1198:   if (!(sp->allNodeValues)) PetscCall(MatCreateVecs(allMat, &(sp->allNodeValues), NULL));
1199:   pointValues = sp->allNodeValues;
1200:   if (!(sp->allDofValues)) PetscCall(MatCreateVecs(allMat, NULL, &(sp->allDofValues)));
1201:   dofValues = sp->allDofValues;
1202:   PetscCall(VecPlaceArray(pointValues, pointEval));
1203:   PetscCall(VecPlaceArray(dofValues, spValue));
1204:   PetscCall(MatMult(allMat, pointValues, dofValues));
1205:   PetscCall(VecResetArray(dofValues));
1206:   PetscCall(VecResetArray(pointValues));
1207:   PetscFunctionReturn(PETSC_SUCCESS);
1208: }

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

1213:   Input Parameters:
1214: + sp        - The `PetscDualSpace` object
1215: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetInteriorData()`

1217:   Output Parameter:
1218: . spValue - The values of interior dual space functionals

1220:   Level: advanced

1222: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1223: @*/
1224: PetscErrorCode PetscDualSpaceApplyInteriorDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1225: {
1226:   Vec pointValues, dofValues;
1227:   Mat intMat;

1229:   PetscFunctionBegin;
1231:   PetscAssertPointer(pointEval, 2);
1232:   PetscAssertPointer(spValue, 3);
1233:   PetscCall(PetscDualSpaceGetInteriorData(sp, NULL, &intMat));
1234:   if (!(sp->intNodeValues)) PetscCall(MatCreateVecs(intMat, &(sp->intNodeValues), NULL));
1235:   pointValues = sp->intNodeValues;
1236:   if (!(sp->intDofValues)) PetscCall(MatCreateVecs(intMat, NULL, &(sp->intDofValues)));
1237:   dofValues = sp->intDofValues;
1238:   PetscCall(VecPlaceArray(pointValues, pointEval));
1239:   PetscCall(VecPlaceArray(dofValues, spValue));
1240:   PetscCall(MatMult(intMat, pointValues, dofValues));
1241:   PetscCall(VecResetArray(dofValues));
1242:   PetscCall(VecResetArray(pointValues));
1243:   PetscFunctionReturn(PETSC_SUCCESS);
1244: }

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

1249:   Input Parameter:
1250: . sp - The dualspace

1252:   Output Parameters:
1253: + allNodes - A `PetscQuadrature` object containing all evaluation nodes
1254: - allMat   - A `Mat` for the node-to-dof transformation

1256:   Level: advanced

1258: .seealso: `PetscQuadrature`, `PetscDualSpace`, `PetscDualSpaceCreate()`, `Mat`
1259: @*/
1260: PetscErrorCode PetscDualSpaceGetAllData(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1261: {
1262:   PetscFunctionBegin;
1264:   if (allNodes) PetscAssertPointer(allNodes, 2);
1265:   if (allMat) PetscAssertPointer(allMat, 3);
1266:   if ((!sp->allNodes || !sp->allMat) && sp->ops->createalldata) {
1267:     PetscQuadrature qpoints;
1268:     Mat             amat;

1270:     PetscUseTypeMethod(sp, createalldata, &qpoints, &amat);
1271:     PetscCall(PetscQuadratureDestroy(&(sp->allNodes)));
1272:     PetscCall(MatDestroy(&(sp->allMat)));
1273:     sp->allNodes = qpoints;
1274:     sp->allMat   = amat;
1275:   }
1276:   if (allNodes) *allNodes = sp->allNodes;
1277:   if (allMat) *allMat = sp->allMat;
1278:   PetscFunctionReturn(PETSC_SUCCESS);
1279: }

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

1284:   Input Parameter:
1285: . sp - The dualspace

1287:   Output Parameters:
1288: + allNodes - A `PetscQuadrature` object containing all evaluation nodes
1289: - allMat   - A `Mat` for the node-to-dof transformation

1291:   Level: advanced

1293: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`, `Mat`, `PetscQuadrature`
1294: @*/
1295: PetscErrorCode PetscDualSpaceCreateAllDataDefault(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1296: {
1297:   PetscInt        spdim;
1298:   PetscInt        numPoints, offset;
1299:   PetscReal      *points;
1300:   PetscInt        f, dim;
1301:   PetscInt        Nc, nrows, ncols;
1302:   PetscInt        maxNumPoints;
1303:   PetscQuadrature q;
1304:   Mat             A;

1306:   PetscFunctionBegin;
1307:   PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
1308:   PetscCall(PetscDualSpaceGetDimension(sp, &spdim));
1309:   if (!spdim) {
1310:     PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allNodes));
1311:     PetscCall(PetscQuadratureSetData(*allNodes, 0, 0, 0, NULL, NULL));
1312:   }
1313:   nrows = spdim;
1314:   PetscCall(PetscDualSpaceGetFunctional(sp, 0, &q));
1315:   PetscCall(PetscQuadratureGetData(q, &dim, NULL, &numPoints, NULL, NULL));
1316:   maxNumPoints = numPoints;
1317:   for (f = 1; f < spdim; f++) {
1318:     PetscInt Np;

1320:     PetscCall(PetscDualSpaceGetFunctional(sp, f, &q));
1321:     PetscCall(PetscQuadratureGetData(q, NULL, NULL, &Np, NULL, NULL));
1322:     numPoints += Np;
1323:     maxNumPoints = PetscMax(maxNumPoints, Np);
1324:   }
1325:   ncols = numPoints * Nc;
1326:   PetscCall(PetscMalloc1(dim * numPoints, &points));
1327:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, nrows, ncols, maxNumPoints * Nc, NULL, &A));
1328:   for (f = 0, offset = 0; f < spdim; f++) {
1329:     const PetscReal *p, *w;
1330:     PetscInt         Np, i;
1331:     PetscInt         fnc;

1333:     PetscCall(PetscDualSpaceGetFunctional(sp, f, &q));
1334:     PetscCall(PetscQuadratureGetData(q, NULL, &fnc, &Np, &p, &w));
1335:     PetscCheck(fnc == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "functional component mismatch");
1336:     for (i = 0; i < Np * dim; i++) points[offset * dim + i] = p[i];
1337:     for (i = 0; i < Np * Nc; i++) PetscCall(MatSetValue(A, f, offset * Nc, w[i], INSERT_VALUES));
1338:     offset += Np;
1339:   }
1340:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1341:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1342:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, allNodes));
1343:   PetscCall(PetscQuadratureSetData(*allNodes, dim, 0, numPoints, points, NULL));
1344:   *allMat = A;
1345:   PetscFunctionReturn(PETSC_SUCCESS);
1346: }

1348: /*@
1349:   PetscDualSpaceGetInteriorData - Get all quadrature points necessary to compute the interior degrees of freedom from
1350:   this space, as well as the matrix that computes the degrees of freedom from the quadrature
1351:   values.

1353:   Input Parameter:
1354: . sp - The dualspace

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

1362:   Level: advanced

1364:   Notes:
1365:   Degrees of freedom are interior degrees of freedom if they belong (by
1366:   `PetscDualSpaceGetSection()`) to interior points in the references, complementary boundary
1367:   degrees of freedom are marked as constrained in the section returned by
1368:   `PetscDualSpaceGetSection()`).

1370: .seealso: `PetscDualSpace`, `PetscQuadrature`, `Mat`, `PetscDualSpaceCreate()`, `PetscDualSpaceGetDimension()`, `PetscDualSpaceGetNumComponents()`, `PetscQuadratureGetData()`
1371: @*/
1372: PetscErrorCode PetscDualSpaceGetInteriorData(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1373: {
1374:   PetscFunctionBegin;
1376:   if (intNodes) PetscAssertPointer(intNodes, 2);
1377:   if (intMat) PetscAssertPointer(intMat, 3);
1378:   if ((!sp->intNodes || !sp->intMat) && sp->ops->createintdata) {
1379:     PetscQuadrature qpoints;
1380:     Mat             imat;

1382:     PetscUseTypeMethod(sp, createintdata, &qpoints, &imat);
1383:     PetscCall(PetscQuadratureDestroy(&(sp->intNodes)));
1384:     PetscCall(MatDestroy(&(sp->intMat)));
1385:     sp->intNodes = qpoints;
1386:     sp->intMat   = imat;
1387:   }
1388:   if (intNodes) *intNodes = sp->intNodes;
1389:   if (intMat) *intMat = sp->intMat;
1390:   PetscFunctionReturn(PETSC_SUCCESS);
1391: }

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

1396:   Input Parameter:
1397: . sp - The dualspace

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

1405:   Level: advanced

1407: .seealso: `PetscDualSpace`, `PetscQuadrature`, `Mat`, `PetscDualSpaceCreate()`, `PetscDualSpaceGetInteriorData()`
1408: @*/
1409: PetscErrorCode PetscDualSpaceCreateInteriorDataDefault(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1410: {
1411:   DM              dm;
1412:   PetscInt        spdim0;
1413:   PetscInt        Nc;
1414:   PetscInt        pStart, pEnd, p, f;
1415:   PetscSection    section;
1416:   PetscInt        numPoints, offset, matoffset;
1417:   PetscReal      *points;
1418:   PetscInt        dim;
1419:   PetscInt       *nnz;
1420:   PetscQuadrature q;
1421:   Mat             imat;

1423:   PetscFunctionBegin;
1425:   PetscCall(PetscDualSpaceGetSection(sp, &section));
1426:   PetscCall(PetscSectionGetConstrainedStorageSize(section, &spdim0));
1427:   if (!spdim0) {
1428:     *intNodes = NULL;
1429:     *intMat   = NULL;
1430:     PetscFunctionReturn(PETSC_SUCCESS);
1431:   }
1432:   PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
1433:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
1434:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
1435:   PetscCall(DMGetDimension(dm, &dim));
1436:   PetscCall(PetscMalloc1(spdim0, &nnz));
1437:   for (p = pStart, f = 0, numPoints = 0; p < pEnd; p++) {
1438:     PetscInt dof, cdof, off, d;

1440:     PetscCall(PetscSectionGetDof(section, p, &dof));
1441:     PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
1442:     if (!(dof - cdof)) continue;
1443:     PetscCall(PetscSectionGetOffset(section, p, &off));
1444:     for (d = 0; d < dof; d++, off++, f++) {
1445:       PetscInt Np;

1447:       PetscCall(PetscDualSpaceGetFunctional(sp, off, &q));
1448:       PetscCall(PetscQuadratureGetData(q, NULL, NULL, &Np, NULL, NULL));
1449:       nnz[f] = Np * Nc;
1450:       numPoints += Np;
1451:     }
1452:   }
1453:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, spdim0, numPoints * Nc, 0, nnz, &imat));
1454:   PetscCall(PetscFree(nnz));
1455:   PetscCall(PetscMalloc1(dim * numPoints, &points));
1456:   for (p = pStart, f = 0, offset = 0, matoffset = 0; p < pEnd; p++) {
1457:     PetscInt dof, cdof, off, d;

1459:     PetscCall(PetscSectionGetDof(section, p, &dof));
1460:     PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
1461:     if (!(dof - cdof)) continue;
1462:     PetscCall(PetscSectionGetOffset(section, p, &off));
1463:     for (d = 0; d < dof; d++, off++, f++) {
1464:       const PetscReal *p;
1465:       const PetscReal *w;
1466:       PetscInt         Np, i;

1468:       PetscCall(PetscDualSpaceGetFunctional(sp, off, &q));
1469:       PetscCall(PetscQuadratureGetData(q, NULL, NULL, &Np, &p, &w));
1470:       for (i = 0; i < Np * dim; i++) points[offset + i] = p[i];
1471:       for (i = 0; i < Np * Nc; i++) PetscCall(MatSetValue(imat, f, matoffset + i, w[i], INSERT_VALUES));
1472:       offset += Np * dim;
1473:       matoffset += Np * Nc;
1474:     }
1475:   }
1476:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, intNodes));
1477:   PetscCall(PetscQuadratureSetData(*intNodes, dim, 0, numPoints, points, NULL));
1478:   PetscCall(MatAssemblyBegin(imat, MAT_FINAL_ASSEMBLY));
1479:   PetscCall(MatAssemblyEnd(imat, MAT_FINAL_ASSEMBLY));
1480:   *intMat = imat;
1481:   PetscFunctionReturn(PETSC_SUCCESS);
1482: }

1484: /*@
1485:   PetscDualSpaceEqual - Determine if two dual spaces are equivalent

1487:   Input Parameters:
1488: + A - A `PetscDualSpace` object
1489: - B - Another `PetscDualSpace` object

1491:   Output Parameter:
1492: . equal - `PETSC_TRUE` if the dual spaces are equivalent

1494:   Level: advanced

1496: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1497: @*/
1498: PetscErrorCode PetscDualSpaceEqual(PetscDualSpace A, PetscDualSpace B, PetscBool *equal)
1499: {
1500:   PetscInt        sizeA, sizeB, dimA, dimB;
1501:   const PetscInt *dofA, *dofB;
1502:   PetscQuadrature quadA, quadB;
1503:   Mat             matA, matB;

1505:   PetscFunctionBegin;
1508:   PetscAssertPointer(equal, 3);
1509:   *equal = PETSC_FALSE;
1510:   PetscCall(PetscDualSpaceGetDimension(A, &sizeA));
1511:   PetscCall(PetscDualSpaceGetDimension(B, &sizeB));
1512:   if (sizeB != sizeA) PetscFunctionReturn(PETSC_SUCCESS);
1513:   PetscCall(DMGetDimension(A->dm, &dimA));
1514:   PetscCall(DMGetDimension(B->dm, &dimB));
1515:   if (dimA != dimB) PetscFunctionReturn(PETSC_SUCCESS);

1517:   PetscCall(PetscDualSpaceGetNumDof(A, &dofA));
1518:   PetscCall(PetscDualSpaceGetNumDof(B, &dofB));
1519:   for (PetscInt d = 0; d < dimA; d++) {
1520:     if (dofA[d] != dofB[d]) PetscFunctionReturn(PETSC_SUCCESS);
1521:   }

1523:   PetscCall(PetscDualSpaceGetInteriorData(A, &quadA, &matA));
1524:   PetscCall(PetscDualSpaceGetInteriorData(B, &quadB, &matB));
1525:   if (!quadA && !quadB) {
1526:     *equal = PETSC_TRUE;
1527:   } else if (quadA && quadB) {
1528:     PetscCall(PetscQuadratureEqual(quadA, quadB, equal));
1529:     if (*equal == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1530:     if (!matA && !matB) PetscFunctionReturn(PETSC_SUCCESS);
1531:     if (matA && matB) PetscCall(MatEqual(matA, matB, equal));
1532:     else *equal = PETSC_FALSE;
1533:   }
1534:   PetscFunctionReturn(PETSC_SUCCESS);
1535: }

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

1540:   Input Parameters:
1541: + sp    - The `PetscDualSpace` object
1542: . f     - The basis functional index
1543: . time  - The time
1544: . cgeom - A context with geometric information for this cell, we currently just use the centroid
1545: . Nc    - The number of components for the function
1546: . func  - The input function
1547: - ctx   - A context for the function

1549:   Output Parameter:
1550: . value - The output value (scalar)

1552:   Calling sequence:
1553: .vb
1554:   PetscErrorCode func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt numComponents, PetscScalar values[], void *ctx)
1555: .ve

1557:   Level: advanced

1559:   Note:
1560:   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.

1562: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1563: @*/
1564: 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)
1565: {
1566:   DM               dm;
1567:   PetscQuadrature  n;
1568:   const PetscReal *points, *weights;
1569:   PetscScalar     *val;
1570:   PetscInt         dimEmbed, qNc, c, Nq, q;

1572:   PetscFunctionBegin;
1574:   PetscAssertPointer(value, 8);
1575:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
1576:   PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
1577:   PetscCall(PetscDualSpaceGetFunctional(sp, f, &n));
1578:   PetscCall(PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights));
1579:   PetscCheck(qNc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_SIZ, "The quadrature components %" PetscInt_FMT " != function components %" PetscInt_FMT, qNc, Nc);
1580:   PetscCall(DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val));
1581:   *value = 0.;
1582:   for (q = 0; q < Nq; ++q) {
1583:     PetscCall((*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx));
1584:     for (c = 0; c < Nc; ++c) *value += val[c] * weights[q * Nc + c];
1585:   }
1586:   PetscCall(DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val));
1587:   PetscFunctionReturn(PETSC_SUCCESS);
1588: }

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

1594:   Not Collective

1596:   Input Parameters:
1597: + sp     - the `PetscDualSpace` object
1598: - height - the height of the mesh point for which the subspace is desired

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

1604:   Level: advanced

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

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

1613: .seealso: `PetscDualSpace`, `PetscSpaceGetHeightSubspace()`, `PetscDualSpaceGetPointSubspace()`
1614: @*/
1615: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1616: {
1617:   PetscInt depth = -1, cStart, cEnd;
1618:   DM       dm;

1620:   PetscFunctionBegin;
1622:   PetscAssertPointer(subsp, 3);
1623:   PetscCheck((sp->uniform), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "A non-uniform dual space does not have a single dual space at each height");
1624:   *subsp = NULL;
1625:   dm     = sp->dm;
1626:   PetscCall(DMPlexGetDepth(dm, &depth));
1627:   PetscCheck(height >= 0 && height <= depth, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid height");
1628:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1629:   if (height == 0 && cEnd == cStart + 1) {
1630:     *subsp = sp;
1631:     PetscFunctionReturn(PETSC_SUCCESS);
1632:   }
1633:   if (!sp->heightSpaces) {
1634:     PetscInt h;
1635:     PetscCall(PetscCalloc1(depth + 1, &(sp->heightSpaces)));

1637:     for (h = 0; h <= depth; h++) {
1638:       if (h == 0 && cEnd == cStart + 1) continue;
1639:       if (sp->ops->createheightsubspace) PetscCall((*sp->ops->createheightsubspace)(sp, height, &(sp->heightSpaces[h])));
1640:       else if (sp->pointSpaces) {
1641:         PetscInt hStart, hEnd;

1643:         PetscCall(DMPlexGetHeightStratum(dm, h, &hStart, &hEnd));
1644:         if (hEnd > hStart) {
1645:           const char *name;

1647:           PetscCall(PetscObjectReference((PetscObject)(sp->pointSpaces[hStart])));
1648:           if (sp->pointSpaces[hStart]) {
1649:             PetscCall(PetscObjectGetName((PetscObject)sp, &name));
1650:             PetscCall(PetscObjectSetName((PetscObject)sp->pointSpaces[hStart], name));
1651:           }
1652:           sp->heightSpaces[h] = sp->pointSpaces[hStart];
1653:         }
1654:       }
1655:     }
1656:   }
1657:   *subsp = sp->heightSpaces[height];
1658:   PetscFunctionReturn(PETSC_SUCCESS);
1659: }

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

1664:   Not Collective

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

1670:   Output Parameters:
1671: . bdsp - the subspace.

1673:   Level: advanced

1675:   Notes:
1676:   The functionals in the subspace are with respect to the intrinsic geometry of the point,
1677:   which will be of lesser dimension if height > 0.

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

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

1685: .seealso: `PetscDualSpace`, `PetscDualSpaceGetHeightSubspace()`
1686: @*/
1687: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1688: {
1689:   PetscInt pStart = 0, pEnd = 0, cStart, cEnd;
1690:   DM       dm;

1692:   PetscFunctionBegin;
1694:   PetscAssertPointer(bdsp, 3);
1695:   *bdsp = NULL;
1696:   dm    = sp->dm;
1697:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1698:   PetscCheck(point >= pStart && point <= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid point");
1699:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1700:   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 */
1701:     *bdsp = sp;
1702:     PetscFunctionReturn(PETSC_SUCCESS);
1703:   }
1704:   if (!sp->pointSpaces) {
1705:     PetscInt p;
1706:     PetscCall(PetscCalloc1(pEnd - pStart, &(sp->pointSpaces)));

1708:     for (p = 0; p < pEnd - pStart; p++) {
1709:       if (p + pStart == cStart && cEnd == cStart + 1) continue;
1710:       if (sp->ops->createpointsubspace) PetscCall((*sp->ops->createpointsubspace)(sp, p + pStart, &(sp->pointSpaces[p])));
1711:       else if (sp->heightSpaces || sp->ops->createheightsubspace) {
1712:         PetscInt dim, depth, height;
1713:         DMLabel  label;

1715:         PetscCall(DMPlexGetDepth(dm, &dim));
1716:         PetscCall(DMPlexGetDepthLabel(dm, &label));
1717:         PetscCall(DMLabelGetValue(label, p + pStart, &depth));
1718:         height = dim - depth;
1719:         PetscCall(PetscDualSpaceGetHeightSubspace(sp, height, &(sp->pointSpaces[p])));
1720:         PetscCall(PetscObjectReference((PetscObject)sp->pointSpaces[p]));
1721:       }
1722:     }
1723:   }
1724:   *bdsp = sp->pointSpaces[point - pStart];
1725:   PetscFunctionReturn(PETSC_SUCCESS);
1726: }

1728: /*@C
1729:   PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis

1731:   Not Collective

1733:   Input Parameter:
1734: . sp - the `PetscDualSpace` object

1736:   Output Parameters:
1737: + perms - Permutations of the interior degrees of freedom, parameterized by the point orientation
1738: - flips - Sign reversal of the interior degrees of freedom, parameterized by the point orientation

1740:   Level: developer

1742:   Note:
1743:   The permutation and flip arrays are organized in the following way
1744: .vb
1745:   perms[p][ornt][dof # on point] = new local dof #
1746:   flips[p][ornt][dof # on point] = reversal or not
1747: .ve

1749: .seealso: `PetscDualSpace`
1750: @*/
1751: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1752: {
1753:   PetscFunctionBegin;
1755:   if (perms) {
1756:     PetscAssertPointer(perms, 2);
1757:     *perms = NULL;
1758:   }
1759:   if (flips) {
1760:     PetscAssertPointer(flips, 3);
1761:     *flips = NULL;
1762:   }
1763:   if (sp->ops->getsymmetries) PetscCall((sp->ops->getsymmetries)(sp, perms, flips));
1764:   PetscFunctionReturn(PETSC_SUCCESS);
1765: }

1767: /*@
1768:   PetscDualSpaceGetFormDegree - Get the form degree k for the k-form the describes the pushforwards/pullbacks of this
1769:   dual space's functionals.

1771:   Input Parameter:
1772: . dsp - The `PetscDualSpace`

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

1782:   Level: developer

1784: .seealso: `PetscDualSpace`, `PetscDTAltV`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1785: @*/
1786: PetscErrorCode PetscDualSpaceGetFormDegree(PetscDualSpace dsp, PetscInt *k)
1787: {
1788:   PetscFunctionBeginHot;
1790:   PetscAssertPointer(k, 2);
1791:   *k = dsp->k;
1792:   PetscFunctionReturn(PETSC_SUCCESS);
1793: }

1795: /*@
1796:   PetscDualSpaceSetFormDegree - Set the form degree k for the k-form the describes the pushforwards/pullbacks of this
1797:   dual space's functionals.

1799:   Input Parameters:
1800: + dsp - The `PetscDualSpace`
1801: - k   - The *signed* degree k of the k.  If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1802:         in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms.  So for example,
1803:         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).
1804:         If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1805:         Hodge star map.  So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1806:         but are stored as 1-forms.

1808:   Level: developer

1810: .seealso: `PetscDualSpace`, `PetscDTAltV`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1811: @*/
1812: PetscErrorCode PetscDualSpaceSetFormDegree(PetscDualSpace dsp, PetscInt k)
1813: {
1814:   PetscInt dim;

1816:   PetscFunctionBeginHot;
1818:   PetscCheck(!dsp->setupcalled, PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up");
1819:   dim = dsp->dm->dim;
1820:   PetscCheck((k >= -dim && k <= dim) || k == PETSC_FORM_DEGREE_UNDEFINED, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported %" PetscInt_FMT "-form on %" PetscInt_FMT "-dimensional reference cell", PetscAbsInt(k), dim);
1821:   dsp->k = k;
1822:   PetscFunctionReturn(PETSC_SUCCESS);
1823: }

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

1828:   Input Parameter:
1829: . dsp - The `PetscDualSpace`

1831:   Output Parameter:
1832: . k - The simplex dimension

1834:   Level: developer

1836:   Note:
1837:   Currently supported values are
1838: .vb
1839:   0: These are H_1 methods that only transform coordinates
1840:   1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1841:   2: These are the same as 1
1842:   3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)
1843: .ve

1845: .seealso: `PetscDualSpace`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1846: @*/
1847: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1848: {
1849:   PetscInt dim;

1851:   PetscFunctionBeginHot;
1853:   PetscAssertPointer(k, 2);
1854:   dim = dsp->dm->dim;
1855:   if (!dsp->k) *k = IDENTITY_TRANSFORM;
1856:   else if (dsp->k == 1) *k = COVARIANT_PIOLA_TRANSFORM;
1857:   else if (dsp->k == -(dim - 1)) *k = CONTRAVARIANT_PIOLA_TRANSFORM;
1858:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported transformation");
1859:   PetscFunctionReturn(PETSC_SUCCESS);
1860: }

1862: /*@C
1863:   PetscDualSpaceTransform - Transform the function values

1865:   Input Parameters:
1866: + dsp       - The `PetscDualSpace`
1867: . trans     - The type of transform
1868: . isInverse - Flag to invert the transform
1869: . fegeom    - The cell geometry
1870: . Nv        - The number of function samples
1871: . Nc        - The number of function components
1872: - vals      - The function values

1874:   Output Parameter:
1875: . vals - The transformed function values

1877:   Level: intermediate

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

1882: .seealso: `PetscDualSpace`, `PetscDualSpaceTransformGradient()`, `PetscDualSpaceTransformHessian()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
1883: @*/
1884: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1885: {
1886:   PetscReal Jstar[9] = {0};
1887:   PetscInt  dim, v, c, Nk;

1889:   PetscFunctionBeginHot;
1891:   PetscAssertPointer(fegeom, 4);
1892:   PetscAssertPointer(vals, 7);
1893:   /* TODO: not handling dimEmbed != dim right now */
1894:   dim = dsp->dm->dim;
1895:   /* No change needed for 0-forms */
1896:   if (!dsp->k) PetscFunctionReturn(PETSC_SUCCESS);
1897:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(dsp->k), &Nk));
1898:   /* TODO: use fegeom->isAffine */
1899:   PetscCall(PetscDTAltVPullbackMatrix(dim, dim, isInverse ? fegeom->J : fegeom->invJ, dsp->k, Jstar));
1900:   for (v = 0; v < Nv; ++v) {
1901:     switch (Nk) {
1902:     case 1:
1903:       for (c = 0; c < Nc; c++) vals[v * Nc + c] *= Jstar[0];
1904:       break;
1905:     case 2:
1906:       for (c = 0; c < Nc; c += 2) DMPlex_Mult2DReal_Internal(Jstar, 1, &vals[v * Nc + c], &vals[v * Nc + c]);
1907:       break;
1908:     case 3:
1909:       for (c = 0; c < Nc; c += 3) DMPlex_Mult3DReal_Internal(Jstar, 1, &vals[v * Nc + c], &vals[v * Nc + c]);
1910:       break;
1911:     default:
1912:       SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported form size %" PetscInt_FMT " for transformation", Nk);
1913:     }
1914:   }
1915:   PetscFunctionReturn(PETSC_SUCCESS);
1916: }

1918: /*@C
1919:   PetscDualSpaceTransformGradient - Transform the function gradient values

1921:   Input Parameters:
1922: + dsp       - The `PetscDualSpace`
1923: . trans     - The type of transform
1924: . isInverse - Flag to invert the transform
1925: . fegeom    - The cell geometry
1926: . Nv        - The number of function gradient samples
1927: . Nc        - The number of function components
1928: - vals      - The function gradient values

1930:   Output Parameter:
1931: . vals - The transformed function gradient values

1933:   Level: intermediate

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

1938: .seealso: `PetscDualSpace`, `PetscDualSpaceTransform()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
1939: @*/
1940: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1941: {
1942:   const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
1943:   PetscInt       v, c, d;

1945:   PetscFunctionBeginHot;
1947:   PetscAssertPointer(fegeom, 4);
1948:   PetscAssertPointer(vals, 7);
1949: #ifdef PETSC_USE_DEBUG
1950:   PetscCheck(dE > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %" PetscInt_FMT, dE);
1951: #endif
1952:   /* Transform gradient */
1953:   if (dim == dE) {
1954:     for (v = 0; v < Nv; ++v) {
1955:       for (c = 0; c < Nc; ++c) {
1956:         switch (dim) {
1957:         case 1:
1958:           vals[(v * Nc + c) * dim] *= fegeom->invJ[0];
1959:           break;
1960:         case 2:
1961:           DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v * Nc + c) * dim], &vals[(v * Nc + c) * dim]);
1962:           break;
1963:         case 3:
1964:           DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v * Nc + c) * dim], &vals[(v * Nc + c) * dim]);
1965:           break;
1966:         default:
1967:           SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1968:         }
1969:       }
1970:     }
1971:   } else {
1972:     for (v = 0; v < Nv; ++v) {
1973:       for (c = 0; c < Nc; ++c) DMPlex_MultTransposeReal_Internal(fegeom->invJ, dim, dE, 1, &vals[(v * Nc + c) * dE], &vals[(v * Nc + c) * dE]);
1974:     }
1975:   }
1976:   /* Assume its a vector, otherwise assume its a bunch of scalars */
1977:   if (Nc == 1 || Nc != dim) PetscFunctionReturn(PETSC_SUCCESS);
1978:   switch (trans) {
1979:   case IDENTITY_TRANSFORM:
1980:     break;
1981:   case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1982:     if (isInverse) {
1983:       for (v = 0; v < Nv; ++v) {
1984:         for (d = 0; d < dim; ++d) {
1985:           switch (dim) {
1986:           case 2:
1987:             DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1988:             break;
1989:           case 3:
1990:             DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1991:             break;
1992:           default:
1993:             SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1994:           }
1995:         }
1996:       }
1997:     } else {
1998:       for (v = 0; v < Nv; ++v) {
1999:         for (d = 0; d < dim; ++d) {
2000:           switch (dim) {
2001:           case 2:
2002:             DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
2003:             break;
2004:           case 3:
2005:             DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
2006:             break;
2007:           default:
2008:             SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
2009:           }
2010:         }
2011:       }
2012:     }
2013:     break;
2014:   case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
2015:     if (isInverse) {
2016:       for (v = 0; v < Nv; ++v) {
2017:         for (d = 0; d < dim; ++d) {
2018:           switch (dim) {
2019:           case 2:
2020:             DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
2021:             break;
2022:           case 3:
2023:             DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
2024:             break;
2025:           default:
2026:             SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
2027:           }
2028:           for (c = 0; c < Nc; ++c) vals[(v * Nc + c) * dim + d] *= fegeom->detJ[0];
2029:         }
2030:       }
2031:     } else {
2032:       for (v = 0; v < Nv; ++v) {
2033:         for (d = 0; d < dim; ++d) {
2034:           switch (dim) {
2035:           case 2:
2036:             DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
2037:             break;
2038:           case 3:
2039:             DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
2040:             break;
2041:           default:
2042:             SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
2043:           }
2044:           for (c = 0; c < Nc; ++c) vals[(v * Nc + c) * dim + d] /= fegeom->detJ[0];
2045:         }
2046:       }
2047:     }
2048:     break;
2049:   }
2050:   PetscFunctionReturn(PETSC_SUCCESS);
2051: }

2053: /*@C
2054:   PetscDualSpaceTransformHessian - Transform the function Hessian values

2056:   Input Parameters:
2057: + dsp       - The `PetscDualSpace`
2058: . trans     - The type of transform
2059: . isInverse - Flag to invert the transform
2060: . fegeom    - The cell geometry
2061: . Nv        - The number of function Hessian samples
2062: . Nc        - The number of function components
2063: - vals      - The function gradient values

2065:   Output Parameter:
2066: . vals - The transformed function Hessian values

2068:   Level: intermediate

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

2073: .seealso: `PetscDualSpace`, `PetscDualSpaceTransform()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
2074: @*/
2075: PetscErrorCode PetscDualSpaceTransformHessian(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
2076: {
2077:   const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
2078:   PetscInt       v, c;

2080:   PetscFunctionBeginHot;
2082:   PetscAssertPointer(fegeom, 4);
2083:   PetscAssertPointer(vals, 7);
2084: #ifdef PETSC_USE_DEBUG
2085:   PetscCheck(dE > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %" PetscInt_FMT, dE);
2086: #endif
2087:   /* Transform Hessian: J^{-T}_{ik} J^{-T}_{jl} H(f)_{kl} = J^{-T}_{ik} H(f)_{kl} J^{-1}_{lj} */
2088:   if (dim == dE) {
2089:     for (v = 0; v < Nv; ++v) {
2090:       for (c = 0; c < Nc; ++c) {
2091:         switch (dim) {
2092:         case 1:
2093:           vals[(v * Nc + c) * dim * dim] *= PetscSqr(fegeom->invJ[0]);
2094:           break;
2095:         case 2:
2096:           DMPlex_PTAP2DReal_Internal(fegeom->invJ, &vals[(v * Nc + c) * dim * dim], &vals[(v * Nc + c) * dim * dim]);
2097:           break;
2098:         case 3:
2099:           DMPlex_PTAP3DReal_Internal(fegeom->invJ, &vals[(v * Nc + c) * dim * dim], &vals[(v * Nc + c) * dim * dim]);
2100:           break;
2101:         default:
2102:           SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
2103:         }
2104:       }
2105:     }
2106:   } else {
2107:     for (v = 0; v < Nv; ++v) {
2108:       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]);
2109:     }
2110:   }
2111:   /* Assume its a vector, otherwise assume its a bunch of scalars */
2112:   if (Nc == 1 || Nc != dim) PetscFunctionReturn(PETSC_SUCCESS);
2113:   switch (trans) {
2114:   case IDENTITY_TRANSFORM:
2115:     break;
2116:   case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
2117:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
2118:   case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
2119:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
2120:   }
2121:   PetscFunctionReturn(PETSC_SUCCESS);
2122: }

2124: /*@C
2125:   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.

2127:   Input Parameters:
2128: + dsp       - The `PetscDualSpace`
2129: . fegeom    - The geometry for this cell
2130: . Nq        - The number of function samples
2131: . Nc        - The number of function components
2132: - pointEval - The function values

2134:   Output Parameter:
2135: . pointEval - The transformed function values

2137:   Level: advanced

2139:   Notes:
2140:   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.

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

2144: .seealso: `PetscDualSpace`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2145: @*/
2146: PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2147: {
2148:   PetscDualSpaceTransformType trans;
2149:   PetscInt                    k;

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

2176: /*@C
2177:   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.

2179:   Input Parameters:
2180: + dsp       - The `PetscDualSpace`
2181: . fegeom    - The geometry for this cell
2182: . Nq        - The number of function samples
2183: . Nc        - The number of function components
2184: - pointEval - The function values

2186:   Output Parameter:
2187: . pointEval - The transformed function values

2189:   Level: advanced

2191:   Notes:
2192:   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.

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

2196: .seealso: `PetscDualSpace`, `PetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2197: @*/
2198: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2199: {
2200:   PetscDualSpaceTransformType trans;
2201:   PetscInt                    k;

2203:   PetscFunctionBeginHot;
2205:   PetscAssertPointer(fegeom, 2);
2206:   PetscAssertPointer(pointEval, 5);
2207:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2208:      This determines their transformation properties. */
2209:   PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
2210:   switch (k) {
2211:   case 0: /* H^1 point evaluations */
2212:     trans = IDENTITY_TRANSFORM;
2213:     break;
2214:   case 1: /* Hcurl preserves tangential edge traces  */
2215:     trans = COVARIANT_PIOLA_TRANSFORM;
2216:     break;
2217:   case 2:
2218:   case 3: /* Hdiv preserve normal traces */
2219:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2220:     break;
2221:   default:
2222:     SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2223:   }
2224:   PetscCall(PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval));
2225:   PetscFunctionReturn(PETSC_SUCCESS);
2226: }

2228: /*@C
2229:   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.

2231:   Input Parameters:
2232: + dsp       - The `PetscDualSpace`
2233: . fegeom    - The geometry for this cell
2234: . Nq        - The number of function gradient samples
2235: . Nc        - The number of function components
2236: - pointEval - The function gradient values

2238:   Output Parameter:
2239: . pointEval - The transformed function gradient values

2241:   Level: advanced

2243:   Notes:
2244:   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.

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

2248: .seealso: `PetscDualSpace`, `PetscDualSpacePushforward()`, `PPetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2249: @*/
2250: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2251: {
2252:   PetscDualSpaceTransformType trans;
2253:   PetscInt                    k;

2255:   PetscFunctionBeginHot;
2257:   PetscAssertPointer(fegeom, 2);
2258:   PetscAssertPointer(pointEval, 5);
2259:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2260:      This determines their transformation properties. */
2261:   PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
2262:   switch (k) {
2263:   case 0: /* H^1 point evaluations */
2264:     trans = IDENTITY_TRANSFORM;
2265:     break;
2266:   case 1: /* Hcurl preserves tangential edge traces  */
2267:     trans = COVARIANT_PIOLA_TRANSFORM;
2268:     break;
2269:   case 2:
2270:   case 3: /* Hdiv preserve normal traces */
2271:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2272:     break;
2273:   default:
2274:     SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2275:   }
2276:   PetscCall(PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval));
2277:   PetscFunctionReturn(PETSC_SUCCESS);
2278: }

2280: /*@C
2281:   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.

2283:   Input Parameters:
2284: + dsp       - The `PetscDualSpace`
2285: . fegeom    - The geometry for this cell
2286: . Nq        - The number of function Hessian samples
2287: . Nc        - The number of function components
2288: - pointEval - The function gradient values

2290:   Output Parameter:
2291: . pointEval - The transformed function Hessian values

2293:   Level: advanced

2295:   Notes:
2296:   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.

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

2300: .seealso: `PetscDualSpace`, `PetscDualSpacePushforward()`, `PPetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2301: @*/
2302: PetscErrorCode PetscDualSpacePushforwardHessian(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2303: {
2304:   PetscDualSpaceTransformType trans;
2305:   PetscInt                    k;

2307:   PetscFunctionBeginHot;
2309:   PetscAssertPointer(fegeom, 2);
2310:   PetscAssertPointer(pointEval, 5);
2311:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2312:      This determines their transformation properties. */
2313:   PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
2314:   switch (k) {
2315:   case 0: /* H^1 point evaluations */
2316:     trans = IDENTITY_TRANSFORM;
2317:     break;
2318:   case 1: /* Hcurl preserves tangential edge traces  */
2319:     trans = COVARIANT_PIOLA_TRANSFORM;
2320:     break;
2321:   case 2:
2322:   case 3: /* Hdiv preserve normal traces */
2323:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2324:     break;
2325:   default:
2326:     SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2327:   }
2328:   PetscCall(PetscDualSpaceTransformHessian(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval));
2329:   PetscFunctionReturn(PETSC_SUCCESS);
2330: }