Actual source code: spacepoly.c

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

      3: static PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscSpace sp, PetscOptionItems *PetscOptionsObject)
      4: {
      5:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;

      7:   PetscFunctionBegin;
      8:   PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace polynomial options");
      9:   PetscCall(PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL));
     10:   PetscOptionsHeadEnd();
     11:   PetscFunctionReturn(PETSC_SUCCESS);
     12: }

     14: static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v)
     15: {
     16:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;

     18:   PetscFunctionBegin;
     19:   PetscCall(PetscViewerASCIIPrintf(v, "%s space of degree %" PetscInt_FMT "\n", poly->tensor ? "Tensor polynomial" : "Polynomial", sp->degree));
     20:   PetscFunctionReturn(PETSC_SUCCESS);
     21: }

     23: static PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
     24: {
     25:   PetscBool iascii;

     27:   PetscFunctionBegin;
     30:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
     31:   if (iascii) PetscCall(PetscSpacePolynomialView_Ascii(sp, viewer));
     32:   PetscFunctionReturn(PETSC_SUCCESS);
     33: }

     35: static PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
     36: {
     37:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;

     39:   PetscFunctionBegin;
     40:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", NULL));
     41:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialSetTensor_C", NULL));
     42:   if (poly->subspaces) {
     43:     PetscInt d;

     45:     for (d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&poly->subspaces[d]));
     46:   }
     47:   PetscCall(PetscFree(poly->subspaces));
     48:   PetscCall(PetscFree(poly));
     49:   PetscFunctionReturn(PETSC_SUCCESS);
     50: }

     52: static PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
     53: {
     54:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;

     56:   PetscFunctionBegin;
     57:   if (poly->setupCalled) PetscFunctionReturn(PETSC_SUCCESS);
     58:   if (sp->Nv <= 1) poly->tensor = PETSC_FALSE;
     59:   if (sp->Nc != 1) {
     60:     PetscInt    Nc     = sp->Nc;
     61:     PetscBool   tensor = poly->tensor;
     62:     PetscInt    Nv     = sp->Nv;
     63:     PetscInt    degree = sp->degree;
     64:     const char *prefix;
     65:     const char *name;
     66:     char        subname[PETSC_MAX_PATH_LEN];
     67:     PetscSpace  subsp;

     69:     PetscCall(PetscSpaceSetType(sp, PETSCSPACESUM));
     70:     PetscCall(PetscSpaceSumSetNumSubspaces(sp, Nc));
     71:     PetscCall(PetscSpaceSumSetInterleave(sp, PETSC_TRUE, PETSC_FALSE));
     72:     PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp));
     73:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix));
     74:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix));
     75:     PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_"));
     76:     if (((PetscObject)sp)->name) {
     77:       PetscCall(PetscObjectGetName((PetscObject)sp, &name));
     78:       PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s sum component", name));
     79:       PetscCall(PetscObjectSetName((PetscObject)subsp, subname));
     80:     } else PetscCall(PetscObjectSetName((PetscObject)subsp, "sum component"));
     81:     PetscCall(PetscSpaceSetType(subsp, PETSCSPACEPOLYNOMIAL));
     82:     PetscCall(PetscSpaceSetDegree(subsp, degree, PETSC_DETERMINE));
     83:     PetscCall(PetscSpaceSetNumComponents(subsp, 1));
     84:     PetscCall(PetscSpaceSetNumVariables(subsp, Nv));
     85:     PetscCall(PetscSpacePolynomialSetTensor(subsp, tensor));
     86:     PetscCall(PetscSpaceSetUp(subsp));
     87:     for (PetscInt i = 0; i < Nc; i++) PetscCall(PetscSpaceSumSetSubspace(sp, i, subsp));
     88:     PetscCall(PetscSpaceDestroy(&subsp));
     89:     PetscCall(PetscSpaceSetUp(sp));
     90:     PetscFunctionReturn(PETSC_SUCCESS);
     91:   }
     92:   if (poly->tensor) {
     93:     sp->maxDegree = PETSC_DETERMINE;
     94:     PetscCall(PetscSpaceSetType(sp, PETSCSPACETENSOR));
     95:     PetscCall(PetscSpaceSetUp(sp));
     96:     PetscFunctionReturn(PETSC_SUCCESS);
     97:   }
     98:   PetscCheck(sp->degree >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %" PetscInt_FMT " invalid", sp->degree);
     99:   sp->maxDegree     = sp->degree;
    100:   poly->setupCalled = PETSC_TRUE;
    101:   PetscFunctionReturn(PETSC_SUCCESS);
    102: }

    104: static PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
    105: {
    106:   PetscInt deg = sp->degree;
    107:   PetscInt n   = sp->Nv;

    109:   PetscFunctionBegin;
    110:   PetscCall(PetscDTBinomialInt(n + deg, n, dim));
    111:   *dim *= sp->Nc;
    112:   PetscFunctionReturn(PETSC_SUCCESS);
    113: }

    115: static PetscErrorCode CoordinateBasis(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt jet, PetscInt Njet, PetscReal pScalar[])
    116: {
    117:   PetscFunctionBegin;
    118:   PetscCall(PetscArrayzero(pScalar, (1 + dim) * Njet * npoints));
    119:   for (PetscInt b = 0; b < 1 + dim; b++) {
    120:     for (PetscInt j = 0; j < PetscMin(1 + dim, Njet); j++) {
    121:       if (j == 0) {
    122:         if (b == 0) {
    123:           for (PetscInt pt = 0; pt < npoints; pt++) pScalar[b * Njet * npoints + j * npoints + pt] = 1.;
    124:         } else {
    125:           for (PetscInt pt = 0; pt < npoints; pt++) pScalar[b * Njet * npoints + j * npoints + pt] = points[pt * dim + (b - 1)];
    126:         }
    127:       } else if (j == b) {
    128:         for (PetscInt pt = 0; pt < npoints; pt++) pScalar[b * Njet * npoints + j * npoints + pt] = 1.;
    129:       }
    130:     }
    131:   }
    132:   PetscFunctionReturn(PETSC_SUCCESS);
    133: }

    135: static PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
    136: {
    137:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
    138:   DM               dm   = sp->dm;
    139:   PetscInt         dim  = sp->Nv;
    140:   PetscInt         Nb, jet, Njet;
    141:   PetscReal       *pScalar;

    143:   PetscFunctionBegin;
    144:   if (!poly->setupCalled) {
    145:     PetscCall(PetscSpaceSetUp(sp));
    146:     PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
    147:     PetscFunctionReturn(PETSC_SUCCESS);
    148:   }
    149:   PetscCheck(!poly->tensor && sp->Nc == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "tensor and multicomponent spaces should have been converted");
    150:   PetscCall(PetscDTBinomialInt(dim + sp->degree, dim, &Nb));
    151:   if (H) {
    152:     jet = 2;
    153:   } else if (D) {
    154:     jet = 1;
    155:   } else {
    156:     jet = 0;
    157:   }
    158:   PetscCall(PetscDTBinomialInt(dim + jet, dim, &Njet));
    159:   PetscCall(DMGetWorkArray(dm, Nb * Njet * npoints, MPIU_REAL, &pScalar));
    160:   // Why are we handling the case degree == 1 specially?  Because we don't want numerical noise when we evaluate hat
    161:   // functions at the vertices of a simplex, which happens when we invert the Vandermonde matrix of the PKD basis.
    162:   // We don't make any promise about which basis is used.
    163:   if (sp->degree == 1) {
    164:     PetscCall(CoordinateBasis(dim, npoints, points, jet, Njet, pScalar));
    165:   } else {
    166:     PetscCall(PetscDTPKDEvalJet(dim, npoints, points, sp->degree, jet, pScalar));
    167:   }
    168:   if (B) {
    169:     PetscInt p_strl = Nb;
    170:     PetscInt b_strl = 1;

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

    175:     PetscCall(PetscArrayzero(B, npoints * Nb));
    176:     for (PetscInt b = 0; b < Nb; b++) {
    177:       for (PetscInt p = 0; p < npoints; p++) B[p * p_strl + b * b_strl] = pScalar[b * b_strr + p * p_strr];
    178:     }
    179:   }
    180:   if (D) {
    181:     PetscInt p_strl = dim * Nb;
    182:     PetscInt b_strl = dim;
    183:     PetscInt d_strl = 1;

    185:     PetscInt b_strr = Njet * npoints;
    186:     PetscInt d_strr = npoints;
    187:     PetscInt p_strr = 1;

    189:     PetscCall(PetscArrayzero(D, npoints * Nb * dim));
    190:     for (PetscInt d = 0; d < dim; d++) {
    191:       for (PetscInt b = 0; b < Nb; b++) {
    192:         for (PetscInt p = 0; p < npoints; p++) D[p * p_strl + b * b_strl + d * d_strl] = pScalar[b * b_strr + (1 + d) * d_strr + p * p_strr];
    193:       }
    194:     }
    195:   }
    196:   if (H) {
    197:     PetscInt p_strl  = dim * dim * Nb;
    198:     PetscInt b_strl  = dim * dim;
    199:     PetscInt d1_strl = dim;
    200:     PetscInt d2_strl = 1;

    202:     PetscInt b_strr = Njet * npoints;
    203:     PetscInt j_strr = npoints;
    204:     PetscInt p_strr = 1;

    206:     PetscInt *derivs;
    207:     PetscCall(PetscCalloc1(dim, &derivs));
    208:     PetscCall(PetscArrayzero(H, npoints * Nb * dim * dim));
    209:     for (PetscInt d1 = 0; d1 < dim; d1++) {
    210:       for (PetscInt d2 = 0; d2 < dim; d2++) {
    211:         PetscInt j;
    212:         derivs[d1]++;
    213:         derivs[d2]++;
    214:         PetscCall(PetscDTGradedOrderToIndex(dim, derivs, &j));
    215:         derivs[d1]--;
    216:         derivs[d2]--;
    217:         for (PetscInt b = 0; b < Nb; b++) {
    218:           for (PetscInt p = 0; p < npoints; p++) H[p * p_strl + b * b_strl + d1 * d1_strl + d2 * d2_strl] = pScalar[b * b_strr + j * j_strr + p * p_strr];
    219:         }
    220:       }
    221:     }
    222:     PetscCall(PetscFree(derivs));
    223:   }
    224:   PetscCall(DMRestoreWorkArray(dm, Nb * Njet * npoints, MPIU_REAL, &pScalar));
    225:   PetscFunctionReturn(PETSC_SUCCESS);
    226: }

    228: /*@
    229:   PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials.

    231:   Input Parameters:
    232: + sp     - the function space object
    233: - tensor - `PETSC_TRUE` for a tensor polynomial space, `PETSC_FALSE` for a polynomial space

    235:   Options Database Key:
    236: . -petscspace_poly_tensor <bool> - Whether to use tensor product polynomials in higher dimension

    238:   Level: intermediate

    240:   Notes:
    241:   It is a tensor space if it is spanned by polynomials whose degree in each variable is
    242:   bounded by the given order, as opposed to the space spanned by polynomials
    243:   whose total degree---summing over all variables---is bounded by the given order.

    245: .seealso: `PetscSpace`, `PetscSpacePolynomialGetTensor()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
    246: @*/
    247: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
    248: {
    249:   PetscFunctionBegin;
    251:   PetscTryMethod(sp, "PetscSpacePolynomialSetTensor_C", (PetscSpace, PetscBool), (sp, tensor));
    252:   PetscFunctionReturn(PETSC_SUCCESS);
    253: }

    255: /*@
    256:   PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor
    257:   polynomials.

    259:   Input Parameter:
    260: . sp - the function space object

    262:   Output Parameter:
    263: . tensor - `PETSC_TRUE` for a tensor polynomial space, `PETSC_FALSE` for a polynomial space

    265:   Level: intermediate

    267:   Notes:
    268:   The space is a tensor space if it is spanned by polynomials whose degree in each variable is
    269:   bounded by the given order, as opposed to the space spanned by polynomials
    270:   whose total degree---summing over all variables---is bounded by the given order.

    272: .seealso: `PetscSpace`, `PetscSpacePolynomialSetTensor()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
    273: @*/
    274: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
    275: {
    276:   PetscFunctionBegin;
    278:   PetscAssertPointer(tensor, 2);
    279:   PetscTryMethod(sp, "PetscSpacePolynomialGetTensor_C", (PetscSpace, PetscBool *), (sp, tensor));
    280:   PetscFunctionReturn(PETSC_SUCCESS);
    281: }

    283: static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
    284: {
    285:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;

    287:   PetscFunctionBegin;
    288:   poly->tensor = tensor;
    289:   PetscFunctionReturn(PETSC_SUCCESS);
    290: }

    292: static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor)
    293: {
    294:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;

    296:   PetscFunctionBegin;
    298:   PetscAssertPointer(tensor, 2);
    299:   *tensor = poly->tensor;
    300:   PetscFunctionReturn(PETSC_SUCCESS);
    301: }

    303: static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp)
    304: {
    305:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
    306:   PetscInt         Nc, dim, order;
    307:   PetscBool        tensor;

    309:   PetscFunctionBegin;
    310:   PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
    311:   PetscCall(PetscSpaceGetNumVariables(sp, &dim));
    312:   PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
    313:   PetscCall(PetscSpacePolynomialGetTensor(sp, &tensor));
    314:   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);
    315:   if (!poly->subspaces) PetscCall(PetscCalloc1(dim, &poly->subspaces));
    316:   if (height <= dim) {
    317:     if (!poly->subspaces[height - 1]) {
    318:       PetscSpace  sub;
    319:       const char *name;

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

    339: static PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
    340: {
    341:   PetscFunctionBegin;
    342:   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_Polynomial;
    343:   sp->ops->setup             = PetscSpaceSetUp_Polynomial;
    344:   sp->ops->view              = PetscSpaceView_Polynomial;
    345:   sp->ops->destroy           = PetscSpaceDestroy_Polynomial;
    346:   sp->ops->getdimension      = PetscSpaceGetDimension_Polynomial;
    347:   sp->ops->evaluate          = PetscSpaceEvaluate_Polynomial;
    348:   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial;
    349:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial));
    350:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial));
    351:   PetscFunctionReturn(PETSC_SUCCESS);
    352: }

    354: /*MC
    355:   PETSCSPACEPOLYNOMIAL = "poly" - A `PetscSpace` object that encapsulates a polynomial space, e.g. P1 is the space of
    356:   linear polynomials. The space is replicated for each component.

    358:   Level: intermediate

    360: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
    361: M*/

    363: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
    364: {
    365:   PetscSpace_Poly *poly;

    367:   PetscFunctionBegin;
    369:   PetscCall(PetscNew(&poly));
    370:   sp->data = poly;

    372:   poly->tensor    = PETSC_FALSE;
    373:   poly->subspaces = NULL;

    375:   PetscCall(PetscSpaceInitialize_Polynomial(sp));
    376:   PetscFunctionReturn(PETSC_SUCCESS);
    377: }