Actual source code: spaceptrimmed.c

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

      3: static PetscErrorCode PetscSpaceSetFromOptions_Ptrimmed(PetscSpace sp, PetscOptionItems PetscOptionsObject)
      4: {
      5:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;

      7:   PetscFunctionBegin;
      8:   PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace polynomial options");
      9:   PetscCall(PetscOptionsInt("-petscspace_ptrimmed_form_degree", "form degree of trimmed space", "PetscSpacePTrimmedSetFormDegree", pt->formDegree, &pt->formDegree, NULL));
     10:   PetscOptionsHeadEnd();
     11:   PetscFunctionReturn(PETSC_SUCCESS);
     12: }

     14: static PetscErrorCode PetscSpacePTrimmedView_Ascii(PetscSpace sp, PetscViewer v)
     15: {
     16:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
     17:   PetscInt             f, tdegree;

     19:   PetscFunctionBegin;
     20:   f       = pt->formDegree;
     21:   tdegree = f == 0 ? sp->degree : sp->degree + 1;
     22:   PetscCall(PetscViewerASCIIPrintf(v, "Trimmed polynomials %" PetscInt_FMT "%s-forms of degree %" PetscInt_FMT " (P-%" PetscInt_FMT "/\\%" PetscInt_FMT ")\n", PetscAbsInt(f), f < 0 ? "*" : "", sp->degree, tdegree, PetscAbsInt(f)));
     23:   PetscFunctionReturn(PETSC_SUCCESS);
     24: }

     26: static PetscErrorCode PetscSpaceView_Ptrimmed(PetscSpace sp, PetscViewer viewer)
     27: {
     28:   PetscBool isascii;

     30:   PetscFunctionBegin;
     33:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
     34:   if (isascii) PetscCall(PetscSpacePTrimmedView_Ascii(sp, viewer));
     35:   PetscFunctionReturn(PETSC_SUCCESS);
     36: }

     38: static PetscErrorCode PetscSpaceDestroy_Ptrimmed(PetscSpace sp)
     39: {
     40:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;

     42:   PetscFunctionBegin;
     43:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePTrimmedGetFormDegree_C", NULL));
     44:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePTrimmedSetFormDegree_C", NULL));
     45:   if (pt->subspaces) {
     46:     PetscInt d;

     48:     for (d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&pt->subspaces[d]));
     49:   }
     50:   PetscCall(PetscFree(pt->subspaces));
     51:   PetscCall(PetscFree(pt));
     52:   PetscFunctionReturn(PETSC_SUCCESS);
     53: }

     55: static PetscErrorCode PetscSpaceSetUp_Ptrimmed(PetscSpace sp)
     56: {
     57:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
     58:   PetscInt             Nf;

     60:   PetscFunctionBegin;
     61:   if (pt->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
     62:   PetscCheck(pt->formDegree >= -sp->Nv && pt->formDegree <= sp->Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Form degree %" PetscInt_FMT " not in valid range [%" PetscInt_FMT ",%" PetscInt_FMT "]", pt->formDegree, sp->Nv, sp->Nv);
     63:   PetscCall(PetscDTBinomialInt(sp->Nv, PetscAbsInt(pt->formDegree), &Nf));
     64:   if (sp->Nc == PETSC_DETERMINE) sp->Nc = Nf;
     65:   PetscCheck(sp->Nc % Nf == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "Number of components %" PetscInt_FMT " is not a multiple of form dimension %" PetscInt_FMT, sp->Nc, Nf);
     66:   if (sp->Nc != Nf) {
     67:     PetscSpace  subsp;
     68:     PetscInt    nCopies = sp->Nc / Nf;
     69:     PetscInt    Nv, deg, maxDeg;
     70:     PetscInt    formDegree = pt->formDegree;
     71:     const char *prefix;
     72:     const char *name;
     73:     char        subname[PETSC_MAX_PATH_LEN];

     75:     PetscCall(PetscSpaceSetType(sp, PETSCSPACESUM));
     76:     PetscCall(PetscSpaceSumSetConcatenate(sp, PETSC_TRUE));
     77:     PetscCall(PetscSpaceSumSetNumSubspaces(sp, nCopies));
     78:     PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp));
     79:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix));
     80:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix));
     81:     PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_"));
     82:     if (((PetscObject)sp)->name) {
     83:       PetscCall(PetscObjectGetName((PetscObject)sp, &name));
     84:       PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s sum component", name));
     85:       PetscCall(PetscObjectSetName((PetscObject)subsp, subname));
     86:     } else PetscCall(PetscObjectSetName((PetscObject)subsp, "sum component"));
     87:     PetscCall(PetscSpaceSetType(subsp, PETSCSPACEPTRIMMED));
     88:     PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
     89:     PetscCall(PetscSpaceSetNumVariables(subsp, Nv));
     90:     PetscCall(PetscSpaceSetNumComponents(subsp, Nf));
     91:     PetscCall(PetscSpaceGetDegree(sp, &deg, &maxDeg));
     92:     PetscCall(PetscSpaceSetDegree(subsp, deg, maxDeg));
     93:     PetscCall(PetscSpacePTrimmedSetFormDegree(subsp, formDegree));
     94:     PetscCall(PetscSpaceSetUp(subsp));
     95:     for (PetscInt i = 0; i < nCopies; i++) PetscCall(PetscSpaceSumSetSubspace(sp, i, subsp));
     96:     PetscCall(PetscSpaceDestroy(&subsp));
     97:     PetscCall(PetscSpaceSetUp(sp));
     98:     PetscFunctionReturn(PETSC_SUCCESS);
     99:   }
    100:   if (sp->degree == PETSC_DEFAULT) sp->degree = 0;
    101:   else PetscCheck(sp->degree >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Invalid negative degree %" PetscInt_FMT, sp->degree);
    102:   sp->maxDegree = (pt->formDegree == 0 || PetscAbsInt(pt->formDegree) == sp->Nv) ? sp->degree : sp->degree + 1;
    103:   if (pt->formDegree == 0 || PetscAbsInt(pt->formDegree) == sp->Nv) {
    104:     // Convert to regular polynomial space
    105:     PetscCall(PetscSpaceSetType(sp, PETSCSPACEPOLYNOMIAL));
    106:     PetscCall(PetscSpaceSetUp(sp));
    107:     PetscFunctionReturn(PETSC_SUCCESS);
    108:   }
    109:   pt->setupcalled = PETSC_TRUE;
    110:   PetscFunctionReturn(PETSC_SUCCESS);
    111: }

    113: static PetscErrorCode PetscSpaceGetDimension_Ptrimmed(PetscSpace sp, PetscInt *dim)
    114: {
    115:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
    116:   PetscInt             f;
    117:   PetscInt             Nf;

    119:   PetscFunctionBegin;
    120:   f = pt->formDegree;
    121:   // For PetscSpace, degree refers to the largest complete polynomial degree contained in the space which
    122:   // is equal to the index of a P trimmed space only for 0-forms: otherwise, the index is degree + 1
    123:   PetscCall(PetscDTPTrimmedSize(sp->Nv, f == 0 ? sp->degree : sp->degree + 1, pt->formDegree, dim));
    124:   PetscCall(PetscDTBinomialInt(sp->Nv, PetscAbsInt(pt->formDegree), &Nf));
    125:   *dim *= (sp->Nc / Nf);
    126:   PetscFunctionReturn(PETSC_SUCCESS);
    127: }

    129: /*
    130:   p in [0, npoints), i in [0, pdim), c in [0, Nc)

    132:   B[p][i][c] = B[p][i_scalar][c][c]
    133: */
    134: static PetscErrorCode PetscSpaceEvaluate_Ptrimmed(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
    135: {
    136:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
    137:   DM                   dm = sp->dm;
    138:   PetscInt             jet, degree, Nf, Ncopies, Njet;
    139:   PetscInt             Nc = sp->Nc;
    140:   PetscInt             f;
    141:   PetscInt             dim = sp->Nv;
    142:   PetscReal           *eval;
    143:   PetscInt             Nb;

    145:   PetscFunctionBegin;
    146:   if (!pt->setupcalled) {
    147:     PetscCall(PetscSpaceSetUp(sp));
    148:     PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
    149:     PetscFunctionReturn(PETSC_SUCCESS);
    150:   }
    151:   if (H) {
    152:     jet = 2;
    153:   } else if (D) {
    154:     jet = 1;
    155:   } else {
    156:     jet = 0;
    157:   }
    158:   f      = pt->formDegree;
    159:   degree = f == 0 ? sp->degree : sp->degree + 1;
    160:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(f), &Nf));
    161:   Ncopies = Nc / Nf;
    162:   PetscCheck(Ncopies == 1, PetscObjectComm((PetscObject)sp), PETSC_ERR_PLIB, "Multicopy spaces should have been converted to PETSCSPACESUM");
    163:   PetscCall(PetscDTBinomialInt(dim + jet, dim, &Njet));
    164:   PetscCall(PetscDTPTrimmedSize(dim, degree, f, &Nb));
    165:   PetscCall(DMGetWorkArray(dm, Nb * Nf * Njet * npoints, MPIU_REAL, &eval));
    166:   PetscCall(PetscDTPTrimmedEvalJet(dim, npoints, points, degree, f, jet, eval));
    167:   if (B) {
    168:     PetscInt p_strl = Nf * Nb;
    169:     PetscInt b_strl = Nf;
    170:     PetscInt v_strl = 1;

    172:     PetscInt b_strr = Nf * Njet * npoints;
    173:     PetscInt v_strr = Njet * npoints;
    174:     PetscInt p_strr = 1;

    176:     for (PetscInt v = 0; v < Nf; v++) {
    177:       for (PetscInt b = 0; b < Nb; b++) {
    178:         for (PetscInt p = 0; p < npoints; p++) B[p * p_strl + b * b_strl + v * v_strl] = eval[b * b_strr + v * v_strr + p * p_strr];
    179:       }
    180:     }
    181:   }
    182:   if (D) {
    183:     PetscInt p_strl = dim * Nf * Nb;
    184:     PetscInt b_strl = dim * Nf;
    185:     PetscInt v_strl = dim;
    186:     PetscInt d_strl = 1;

    188:     PetscInt b_strr = Nf * Njet * npoints;
    189:     PetscInt v_strr = Njet * npoints;
    190:     PetscInt d_strr = npoints;
    191:     PetscInt p_strr = 1;

    193:     for (PetscInt v = 0; v < Nf; v++) {
    194:       for (PetscInt d = 0; d < dim; d++) {
    195:         for (PetscInt b = 0; b < Nb; b++) {
    196:           for (PetscInt p = 0; p < npoints; p++) D[p * p_strl + b * b_strl + v * v_strl + d * d_strl] = eval[b * b_strr + v * v_strr + (1 + d) * d_strr + p * p_strr];
    197:         }
    198:       }
    199:     }
    200:   }
    201:   if (H) {
    202:     PetscInt p_strl  = dim * dim * Nf * Nb;
    203:     PetscInt b_strl  = dim * dim * Nf;
    204:     PetscInt v_strl  = dim * dim;
    205:     PetscInt d1_strl = dim;
    206:     PetscInt d2_strl = 1;

    208:     PetscInt b_strr = Nf * Njet * npoints;
    209:     PetscInt v_strr = Njet * npoints;
    210:     PetscInt j_strr = npoints;
    211:     PetscInt p_strr = 1;

    213:     PetscInt *derivs;
    214:     PetscCall(PetscCalloc1(dim, &derivs));
    215:     for (PetscInt d1 = 0; d1 < dim; d1++) {
    216:       for (PetscInt d2 = 0; d2 < dim; d2++) {
    217:         PetscInt j;
    218:         derivs[d1]++;
    219:         derivs[d2]++;
    220:         PetscCall(PetscDTGradedOrderToIndex(dim, derivs, &j));
    221:         derivs[d1]--;
    222:         derivs[d2]--;
    223:         for (PetscInt v = 0; v < Nf; v++) {
    224:           for (PetscInt b = 0; b < Nb; b++) {
    225:             for (PetscInt p = 0; p < npoints; p++) H[p * p_strl + b * b_strl + v * v_strl + d1 * d1_strl + d2 * d2_strl] = eval[b * b_strr + v * v_strr + j * j_strr + p * p_strr];
    226:           }
    227:         }
    228:       }
    229:     }
    230:     PetscCall(PetscFree(derivs));
    231:   }
    232:   PetscCall(DMRestoreWorkArray(dm, Nb * Nf * Njet * npoints, MPIU_REAL, &eval));
    233:   PetscFunctionReturn(PETSC_SUCCESS);
    234: }

    236: /*@
    237:   PetscSpacePTrimmedSetFormDegree - Set the form degree of the trimmed polynomials.

    239:   Input Parameters:
    240: + sp         - the function space object
    241: - formDegree - the form degree

    243:   Options Database Key:
    244: . -petscspace_ptrimmed_form_degree <int> - The trimmed polynomial form degree

    246:   Level: intermediate

    248: .seealso: `PetscSpace`, `PetscDTAltV`, `PetscDTPTrimmedEvalJet()`, `PetscSpacePTrimmedGetFormDegree()`
    249: @*/
    250: PetscErrorCode PetscSpacePTrimmedSetFormDegree(PetscSpace sp, PetscInt formDegree)
    251: {
    252:   PetscFunctionBegin;
    254:   PetscTryMethod(sp, "PetscSpacePTrimmedSetFormDegree_C", (PetscSpace, PetscInt), (sp, formDegree));
    255:   PetscFunctionReturn(PETSC_SUCCESS);
    256: }

    258: /*@
    259:   PetscSpacePTrimmedGetFormDegree - Get the form degree of the trimmed polynomials.

    261:   Input Parameter:
    262: . sp - the function space object

    264:   Output Parameter:
    265: . formDegree - the form degree

    267:   Level: intermediate

    269: .seealso: `PetscSpace`, `PetscDTAltV`, `PetscDTPTrimmedEvalJet()`, `PetscSpacePTrimmedSetFormDegree()`
    270: @*/
    271: PetscErrorCode PetscSpacePTrimmedGetFormDegree(PetscSpace sp, PetscInt *formDegree)
    272: {
    273:   PetscFunctionBegin;
    275:   PetscAssertPointer(formDegree, 2);
    276:   PetscTryMethod(sp, "PetscSpacePTrimmedGetFormDegree_C", (PetscSpace, PetscInt *), (sp, formDegree));
    277:   PetscFunctionReturn(PETSC_SUCCESS);
    278: }

    280: static PetscErrorCode PetscSpacePTrimmedSetFormDegree_Ptrimmed(PetscSpace sp, PetscInt formDegree)
    281: {
    282:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;

    284:   PetscFunctionBegin;
    285:   pt->formDegree = formDegree;
    286:   PetscFunctionReturn(PETSC_SUCCESS);
    287: }

    289: static PetscErrorCode PetscSpacePTrimmedGetFormDegree_Ptrimmed(PetscSpace sp, PetscInt *formDegree)
    290: {
    291:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;

    293:   PetscFunctionBegin;
    295:   PetscAssertPointer(formDegree, 2);
    296:   *formDegree = pt->formDegree;
    297:   PetscFunctionReturn(PETSC_SUCCESS);
    298: }

    300: static PetscErrorCode PetscSpaceGetHeightSubspace_Ptrimmed(PetscSpace sp, PetscInt height, PetscSpace *subsp)
    301: {
    302:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
    303:   PetscInt             dim;

    305:   PetscFunctionBegin;
    306:   PetscCall(PetscSpaceGetNumVariables(sp, &dim));
    307:   PetscCheck(height <= dim && height >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %" PetscInt_FMT " for dimension %" PetscInt_FMT " space", height, dim);
    308:   if (!pt->subspaces) PetscCall(PetscCalloc1(dim, &pt->subspaces));
    309:   if ((dim - height) <= PetscAbsInt(pt->formDegree)) {
    310:     if (!pt->subspaces[height - 1]) {
    311:       PetscInt    Nc, degree, Nf, Ncopies, Nfsub;
    312:       PetscSpace  sub;
    313:       const char *name;

    315:       PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
    316:       PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(pt->formDegree), &Nf));
    317:       PetscCall(PetscDTBinomialInt(dim - height, PetscAbsInt(pt->formDegree), &Nfsub));
    318:       Ncopies = Nf / Nc;
    319:       PetscCall(PetscSpaceGetDegree(sp, &degree, NULL));

    321:       PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub));
    322:       PetscCall(PetscObjectGetName((PetscObject)sp, &name));
    323:       PetscCall(PetscObjectSetName((PetscObject)sub, name));
    324:       PetscCall(PetscSpaceSetType(sub, PETSCSPACEPTRIMMED));
    325:       PetscCall(PetscSpaceSetNumComponents(sub, Nfsub * Ncopies));
    326:       PetscCall(PetscSpaceSetDegree(sub, degree, PETSC_DETERMINE));
    327:       PetscCall(PetscSpaceSetNumVariables(sub, dim - height));
    328:       PetscCall(PetscSpacePTrimmedSetFormDegree(sub, pt->formDegree));
    329:       PetscCall(PetscSpaceSetUp(sub));
    330:       pt->subspaces[height - 1] = sub;
    331:     }
    332:     *subsp = pt->subspaces[height - 1];
    333:   } else {
    334:     *subsp = NULL;
    335:   }
    336:   PetscFunctionReturn(PETSC_SUCCESS);
    337: }

    339: static PetscErrorCode PetscSpaceInitialize_Ptrimmed(PetscSpace sp)
    340: {
    341:   PetscFunctionBegin;
    342:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePTrimmedGetFormDegree_C", PetscSpacePTrimmedGetFormDegree_Ptrimmed));
    343:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePTrimmedSetFormDegree_C", PetscSpacePTrimmedSetFormDegree_Ptrimmed));
    344:   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_Ptrimmed;
    345:   sp->ops->setup             = PetscSpaceSetUp_Ptrimmed;
    346:   sp->ops->view              = PetscSpaceView_Ptrimmed;
    347:   sp->ops->destroy           = PetscSpaceDestroy_Ptrimmed;
    348:   sp->ops->getdimension      = PetscSpaceGetDimension_Ptrimmed;
    349:   sp->ops->evaluate          = PetscSpaceEvaluate_Ptrimmed;
    350:   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Ptrimmed;
    351:   PetscFunctionReturn(PETSC_SUCCESS);
    352: }

    354: /*MC
    355:   PETSCSPACEPTRIMMED = "ptrimmed" - A `PetscSpace` object that encapsulates a trimmed polynomial space.

    357:   Trimmed polynomial spaces are defined for $k$-forms, and are defined by
    358:   $
    359:   \mathcal{P}^-_r \Lambda^k(\mathbb{R}^n) = mathcal{P}_{r-1} \Lambda^k(\mathbb{R}^n) \oplus \kappa [\mathcal{H}_{r-1} \Lambda^{k+1}(\mathbb{R}^n)],
    360:   $
    361:   where $\mathcal{H}_{r-1}$ are homogeneous polynomials and $\kappa$ is the Koszul differential.  This decomposition is detailed in ``Finite element exterior calculus'', Arnold, 2018.

    363:   Level: intermediate

    365:   Notes:
    366:   Trimmed polynomial spaces correspond to several common conformal approximation spaces in the de Rham complex:

    368:   In $H^1$ ($\sim k=0$), trimmed polynomial spaces are identical to the standard polynomial spaces, $\mathcal{P}_r^- \sim P_r$.

    370:   In $H(\text{curl})$, ($\sim k=1$), trimmed polynomial spaces are equivalent to $H(\text{curl})$-Nedelec spaces of the first kind and can be written as
    371:   $
    372:     \begin{cases}
    373:       [P_{r-1}(\mathbb{R}^2)]^2 \oplus \mathrm{rot}(\bf{x}) H_{r-1}(\mathbb{R}^2), & n = 2, \\
    374:       [P_{r-1}(\mathbb{R}^3)]^3 \oplus \bf{x} \times [H_{r-1}(\mathbb{R}^3)]^3, & n = 3.
    375:     \end{cases}
    376:   $

    378:   In $H(\text{div})$ ($\sim k=n-1$), trimmed polynomial spaces are equivalent to Raviart-Thomas spaces ($n=2$) and $H(\text{div})$-Nedelec spaces of the first kind ($n=3$), and can be written as
    379:   $
    380:     [P_{r-1}(\mathbb{R}^n)]^n \oplus \bf{x} H_{r-1}(\mathbb{R}^n).
    381:   $

    383:   In $L_2$, ($\sim k=n$), trimmed polynomial spaces are identical to the standard polynomial spaces of one degree less, $\mathcal{P}_r^- \sim P_{r-1}$.

    385: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscDTPTrimmedEvalJet()`
    386: M*/

    388: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Ptrimmed(PetscSpace sp)
    389: {
    390:   PetscSpace_Ptrimmed *pt;

    392:   PetscFunctionBegin;
    394:   PetscCall(PetscNew(&pt));
    395:   sp->data = pt;

    397:   pt->subspaces = NULL;
    398:   sp->Nc        = PETSC_DETERMINE;

    400:   PetscCall(PetscSpaceInitialize_Ptrimmed(sp));
    401:   PetscFunctionReturn(PETSC_SUCCESS);
    402: }