Actual source code: dspacesimple.c

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

      4: static PetscErrorCode PetscDualSpaceSetUp_Simple(PetscDualSpace sp)
      5: {
      6:   PetscDualSpace_Simple *s  = (PetscDualSpace_Simple *)sp->data;
      7:   DM                     dm = sp->dm;
      8:   PetscInt               dim, pStart, pEnd;
      9:   PetscSection           section;

     11:   PetscFunctionBegin;
     12:   PetscCall(DMGetDimension(dm, &dim));
     13:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
     14:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &section));
     15:   PetscCall(PetscSectionSetChart(section, pStart, pEnd));
     16:   PetscCall(PetscSectionSetDof(section, pStart, s->dim));
     17:   PetscCall(PetscSectionSetUp(section));
     18:   sp->pointSection = section;
     19:   PetscFunctionReturn(PETSC_SUCCESS);
     20: }

     22: static PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
     23: {
     24:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;

     26:   PetscFunctionBegin;
     27:   PetscCall(PetscFree(s->numDof));
     28:   PetscCall(PetscFree(s));
     29:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetDimension_C", NULL));
     30:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetFunctional_C", NULL));
     31:   PetscFunctionReturn(PETSC_SUCCESS);
     32: }

     34: static PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace spNew)
     35: {
     36:   PetscInt dim, d;

     38:   PetscFunctionBegin;
     39:   PetscCall(PetscDualSpaceGetDimension(sp, &dim));
     40:   PetscCall(PetscDualSpaceSimpleSetDimension(spNew, dim));
     41:   for (d = 0; d < dim; ++d) {
     42:     PetscQuadrature q;

     44:     PetscCall(PetscDualSpaceGetFunctional(sp, d, &q));
     45:     PetscCall(PetscDualSpaceSimpleSetFunctional(spNew, d, q));
     46:   }
     47:   PetscFunctionReturn(PETSC_SUCCESS);
     48: }

     50: static PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
     51: {
     52:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
     53:   DM                     dm;
     54:   PetscInt               spatialDim, f;

     56:   PetscFunctionBegin;
     57:   for (f = 0; f < s->dim; ++f) PetscCall(PetscQuadratureDestroy(&sp->functional[f]));
     58:   PetscCall(PetscFree(sp->functional));
     59:   s->dim = dim;
     60:   PetscCall(PetscCalloc1(s->dim, &sp->functional));
     61:   PetscCall(PetscFree(s->numDof));
     62:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
     63:   PetscCall(DMGetCoordinateDim(dm, &spatialDim));
     64:   PetscCall(PetscCalloc1(spatialDim + 1, &s->numDof));
     65:   s->numDof[spatialDim] = dim;
     66:   PetscFunctionReturn(PETSC_SUCCESS);
     67: }

     69: static PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
     70: {
     71:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
     72:   PetscReal             *weights;
     73:   PetscInt               Nc, c, Nq, p;

     75:   PetscFunctionBegin;
     76:   PetscCheck(!(f < 0) && !(f >= s->dim), PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", f, s->dim);
     77:   PetscCall(PetscQuadratureDuplicate(q, &sp->functional[f]));
     78:   /* Reweight so that it has unit volume: Do we want to do this for Nc > 1? */
     79:   PetscCall(PetscQuadratureGetData(sp->functional[f], NULL, &Nc, &Nq, NULL, (const PetscReal **)&weights));
     80:   for (c = 0; c < Nc; ++c) {
     81:     PetscReal vol = 0.0;

     83:     for (p = 0; p < Nq; ++p) vol += weights[p * Nc + c];
     84:     for (p = 0; p < Nq; ++p) weights[p * Nc + c] /= (vol == 0.0 ? 1.0 : vol);
     85:   }
     86:   PetscFunctionReturn(PETSC_SUCCESS);
     87: }

     89: /*@
     90:   PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis

     92:   Logically Collective

     94:   Input Parameters:
     95: + sp  - the `PetscDualSpace`
     96: - dim - the basis dimension

     98:   Level: intermediate

    100: .seealso: `PETSCDUALSPACESIMPLE`, `PetscDualSpace`, `PetscDualSpaceSimpleSetFunctional()`
    101: @*/
    102: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
    103: {
    104:   PetscFunctionBegin;
    107:   PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change dimension after dual space is set up");
    108:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace, PetscInt), (sp, dim));
    109:   PetscFunctionReturn(PETSC_SUCCESS);
    110: }

    112: /*@
    113:   PetscDualSpaceSimpleSetFunctional - Set the given basis functional for this dual space

    115:   Not Collective

    117:   Input Parameters:
    118: + sp   - the `PetscDualSpace`
    119: . func - the basis index
    120: - q    - the basis functional

    122:   Level: intermediate

    124:   Note:
    125:   The quadrature will be reweighted so that it has unit volume.

    127: .seealso: `PETSCDUALSPACESIMPLE`, `PetscDualSpace`, `PetscDualSpaceSimpleSetDimension()`
    128: @*/
    129: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
    130: {
    131:   PetscFunctionBegin;
    133:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace, PetscInt, PetscQuadrature), (sp, func, q));
    134:   PetscFunctionReturn(PETSC_SUCCESS);
    135: }

    137: static PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
    138: {
    139:   PetscFunctionBegin;
    140:   sp->ops->setup                = PetscDualSpaceSetUp_Simple;
    141:   sp->ops->view                 = NULL;
    142:   sp->ops->destroy              = PetscDualSpaceDestroy_Simple;
    143:   sp->ops->duplicate            = PetscDualSpaceDuplicate_Simple;
    144:   sp->ops->createheightsubspace = NULL;
    145:   sp->ops->createpointsubspace  = NULL;
    146:   sp->ops->getsymmetries        = NULL;
    147:   sp->ops->apply                = PetscDualSpaceApplyDefault;
    148:   sp->ops->applyall             = PetscDualSpaceApplyAllDefault;
    149:   sp->ops->applyint             = PetscDualSpaceApplyInteriorDefault;
    150:   sp->ops->createalldata        = PetscDualSpaceCreateAllDataDefault;
    151:   sp->ops->createintdata        = PetscDualSpaceCreateInteriorDataDefault;
    152:   PetscFunctionReturn(PETSC_SUCCESS);
    153: }

    155: /*MC
    156:   PETSCDUALSPACESIMPLE = "simple" - A `PetscDualSpaceType` that encapsulates a dual space of functionals provided with `PetscDualSpaceSimpleSetFunctional()`

    158:   Level: intermediate

    160:   Developer Note:
    161:   It is not clear this has a good name

    163: .seealso: `PetscDualSpace`, `PetscDualSpaceSimpleSetFunctional()`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
    164: M*/

    166: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
    167: {
    168:   PetscDualSpace_Simple *s;

    170:   PetscFunctionBegin;
    172:   PetscCall(PetscNew(&s));
    173:   sp->data = s;

    175:   s->dim    = 0;
    176:   s->numDof = NULL;

    178:   PetscCall(PetscDualSpaceInitialize_Simple(sp));
    179:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple));
    180:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple));
    181:   PetscFunctionReturn(PETSC_SUCCESS);
    182: }