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, No Fortran Support

     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, PeOp PetscObject obj, const char name[])
    224: {
    225:   PetscFunctionBegin;
    227:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
    228:   PetscFunctionReturn(PETSC_SUCCESS);
    229: }

    231: /*@C
    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 isascii;

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

    257: /*@C
    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: /*@C
    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: /*@C
    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: /*@C
    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:   PetscCall(PetscFEInitializePackage());

    452:   PetscCall(PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView));
    453:   s->order       = 0;
    454:   s->Nc          = 1;
    455:   s->k           = 0;
    456:   s->spdim       = -1;
    457:   s->spintdim    = -1;
    458:   s->uniform     = PETSC_TRUE;
    459:   s->setupcalled = PETSC_FALSE;
    460:   *sp            = s;
    461:   PetscFunctionReturn(PETSC_SUCCESS);
    462: }

    464: /*@C
    465:   PetscDualSpaceDuplicate - Creates a duplicate `PetscDualSpace` object that is not setup.

    467:   Collective

    469:   Input Parameter:
    470: . sp - The original `PetscDualSpace`

    472:   Output Parameter:
    473: . spNew - The duplicate `PetscDualSpace`

    475:   Level: beginner

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

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

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

    504: /*@C
    505:   PetscDualSpaceGetDM - Get the `DM` representing the reference cell of a `PetscDualSpace`

    507:   Not Collective

    509:   Input Parameter:
    510: . sp - The `PetscDualSpace`

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

    515:   Level: intermediate

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

    528: /*@C
    529:   PetscDualSpaceSetDM - Get the `DM` representing the reference cell

    531:   Not Collective

    533:   Input Parameters:
    534: + sp - The `PetscDual`Space
    535: - dm - The reference cell

    537:   Level: intermediate

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

    554: /*@C
    555:   PetscDualSpaceGetOrder - Get the order of the dual space

    557:   Not Collective

    559:   Input Parameter:
    560: . sp - The `PetscDualSpace`

    562:   Output Parameter:
    563: . order - The order

    565:   Level: intermediate

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

    578: /*@C
    579:   PetscDualSpaceSetOrder - Set the order of the dual space

    581:   Not Collective

    583:   Input Parameters:
    584: + sp    - The `PetscDualSpace`
    585: - order - The order

    587:   Level: intermediate

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

    600: /*@C
    601:   PetscDualSpaceGetNumComponents - Return the number of components for this space

    603:   Input Parameter:
    604: . sp - The `PetscDualSpace`

    606:   Output Parameter:
    607: . Nc - The number of components

    609:   Level: intermediate

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

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

    625: /*@C
    626:   PetscDualSpaceSetNumComponents - Set the number of components for this space

    628:   Input Parameters:
    629: + sp - The `PetscDualSpace`
    630: - Nc - The number of components

    632:   Level: intermediate

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

    645: /*@C
    646:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

    648:   Not Collective

    650:   Input Parameters:
    651: + sp - The `PetscDualSpace`
    652: - i  - The basis number

    654:   Output Parameter:
    655: . functional - The basis functional

    657:   Level: intermediate

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

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

    674: /*@C
    675:   PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals

    677:   Not Collective

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

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

    685:   Level: intermediate

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

    697:     PetscCall(PetscDualSpaceGetSection(sp, &section));
    698:     if (section) PetscCall(PetscSectionGetStorageSize(section, &sp->spdim));
    699:     else sp->spdim = 0;
    700:   }
    701:   *dim = sp->spdim;
    702:   PetscFunctionReturn(PETSC_SUCCESS);
    703: }

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

    708:   Not Collective

    710:   Input Parameter:
    711: . sp - The `PetscDualSpace`

    713:   Output Parameter:
    714: . intdim - The dimension

    716:   Level: intermediate

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

    728:     PetscCall(PetscDualSpaceGetSection(sp, &section));
    729:     if (section) PetscCall(PetscSectionGetConstrainedStorageSize(section, &sp->spintdim));
    730:     else sp->spintdim = 0;
    731:   }
    732:   *intdim = sp->spintdim;
    733:   PetscFunctionReturn(PETSC_SUCCESS);
    734: }

    736: /*@C
    737:   PetscDualSpaceGetUniform - Whether this dual space is uniform

    739:   Not Collective

    741:   Input Parameter:
    742: . sp - A dual space

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

    748:   Level: advanced

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

    754: .seealso: `PetscDualSpace`, `PetscDualSpaceGetPointSubspace()`, `PetscDualSpaceGetSymmetries()`
    755: @*/
    756: PetscErrorCode PetscDualSpaceGetUniform(PetscDualSpace sp, PetscBool *uniform)
    757: {
    758:   PetscFunctionBegin;
    760:   PetscAssertPointer(uniform, 2);
    761:   *uniform = sp->uniform;
    762:   PetscFunctionReturn(PETSC_SUCCESS);
    763: }

    765: /*@CC
    766:   PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension

    768:   Not Collective

    770:   Input Parameter:
    771: . sp - The `PetscDualSpace`

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

    776:   Level: intermediate

    778:   Note:
    779:   Do not free `numDof`

    781: .seealso: `PetscDualSpace`, `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
    782: @*/
    783: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt *numDof[])
    784: {
    785:   PetscFunctionBegin;
    787:   PetscAssertPointer(numDof, 2);
    788:   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");
    789:   if (!sp->numDof) {
    790:     DM           dm;
    791:     PetscInt     depth, d;
    792:     PetscSection section;

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

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

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

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

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

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

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

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

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

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

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

    892: /*@C
    893:   PetscDualSpaceGetSection - Create a `PetscSection` over the reference cell with the layout from this space

    895:   Collective

    897:   Input Parameter:
    898: . sp - The `PetscDualSpace`

    900:   Output Parameter:
    901: . section - The section

    903:   Level: advanced

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

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

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

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

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

    941:   Collective

    943:   Input Parameter:
    944: . sp - The `PetscDualSpace`

    946:   Output Parameter:
    947: . section - The interior section

    949:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

    1044:   Output Parameter:
    1045: . value - numComp output values

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

    1052:   Level: beginner

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

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

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

    1073:   Output Parameter:
    1074: . spValue - The values of all dual space functionals

    1076:   Level: advanced

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

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

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

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

    1098:   Level: advanced

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

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

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

    1122:   Output Parameter:
    1123: . value - The output value

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

    1130:   Level: advanced

    1132:   Note:
    1133:   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.

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

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

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

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

    1179:   Output Parameter:
    1180: . spValue - The values of all dual space functionals

    1182:   Level: advanced

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

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

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

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

    1215:   Output Parameter:
    1216: . spValue - The values of interior dual space functionals

    1218:   Level: advanced

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

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

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

    1247:   Input Parameter:
    1248: . sp - The dualspace

    1250:   Output Parameters:
    1251: + allNodes - A `PetscQuadrature` object containing all evaluation nodes, pass `NULL` if not needed
    1252: - allMat   - A `Mat` for the node-to-dof transformation, pass `NULL` if not needed

    1254:   Level: advanced

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

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

    1279: /*@C
    1280:   PetscDualSpaceCreateAllDataDefault - Create all evaluation nodes and the node-to-dof matrix by examining functionals

    1282:   Input Parameter:
    1283: . sp - The dualspace

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

    1289:   Level: advanced

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

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

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

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

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

    1351:   Input Parameter:
    1352: . sp - The dualspace

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

    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, PeOp PetscQuadrature *intNodes, PeOp 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: /*@C
    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: /*@C
    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: /*@C
    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) PetscUseTypeMethod(sp, createheightsubspace, 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: /*@C
    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 Parameter:
    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) PetscUseTypeMethod(sp, createpointsubspace, 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:   PetscTryTypeMethod(sp, getsymmetries, perms, flips);
    1764:   PetscFunctionReturn(PETSC_SUCCESS);
    1765: }

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

    2051: /*@C
    2052:   PetscDualSpaceTransformHessian - Transform the function Hessian values

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

    2063:   Output Parameter:
    2064: . vals - The transformed function Hessian values

    2066:   Level: intermediate

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

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

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

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

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

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

    2133:   Level: advanced

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

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

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

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

    2172: /*@C
    2173:   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.

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

    2182:   Output Parameter:
    2183: . pointEval - The transformed function values

    2185:   Level: advanced

    2187:   Notes:
    2188:   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.

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

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

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

    2224: /*@C
    2225:   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.

    2227:   Input Parameters:
    2228: + dsp       - The `PetscDualSpace`
    2229: . fegeom    - The geometry for this cell
    2230: . Nq        - The number of function gradient samples
    2231: . Nc        - The number of function components
    2232: - pointEval - The function gradient values

    2234:   Output Parameter:
    2235: . pointEval - The transformed function gradient values

    2237:   Level: advanced

    2239:   Notes:
    2240:   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.

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

    2244: .seealso: `PetscDualSpace`, `PetscDualSpacePushforward()`, `PetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
    2245: @*/
    2246: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
    2247: {
    2248:   PetscDualSpaceTransformType trans;
    2249:   PetscInt                    k;

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

    2276: /*@C
    2277:   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.

    2279:   Input Parameters:
    2280: + dsp       - The `PetscDualSpace`
    2281: . fegeom    - The geometry for this cell
    2282: . Nq        - The number of function Hessian samples
    2283: . Nc        - The number of function components
    2284: - pointEval - The function gradient values

    2286:   Output Parameter:
    2287: . pointEval - The transformed function Hessian values

    2289:   Level: advanced

    2291:   Notes:
    2292:   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.

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

    2296: .seealso: `PetscDualSpace`, `PetscDualSpacePushforward()`, `PetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
    2297: @*/
    2298: PetscErrorCode PetscDualSpacePushforwardHessian(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
    2299: {
    2300:   PetscDualSpaceTransformType trans;
    2301:   PetscInt                    k;

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