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

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

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

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

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

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

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

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

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

188:     PetscCall(PetscArrayzero(D, npoints * Nb * dim));
189:     for (PetscInt d = 0; d < dim; d++) {
190:       for (PetscInt b = 0; b < Nb; b++) {
191:         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];
192:       }
193:     }
194:   }
195:   if (H) {
196:     PetscInt p_strl  = dim * dim * Nb;
197:     PetscInt b_strl  = dim * dim;
198:     PetscInt d1_strl = dim;
199:     PetscInt d2_strl = 1;

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

205:     PetscInt *derivs;
206:     PetscCall(PetscCalloc1(dim, &derivs));
207:     PetscCall(PetscArrayzero(H, npoints * Nb * dim * dim));
208:     for (PetscInt d1 = 0; d1 < dim; d1++) {
209:       for (PetscInt d2 = 0; d2 < dim; d2++) {
210:         PetscInt j;
211:         derivs[d1]++;
212:         derivs[d2]++;
213:         PetscCall(PetscDTGradedOrderToIndex(dim, derivs, &j));
214:         derivs[d1]--;
215:         derivs[d2]--;
216:         for (PetscInt b = 0; b < Nb; b++) {
217:           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];
218:         }
219:       }
220:     }
221:     PetscCall(PetscFree(derivs));
222:   }
223:   PetscCall(DMRestoreWorkArray(dm, Nb * Njet * npoints, MPIU_REAL, &pScalar));
224:   PetscFunctionReturn(PETSC_SUCCESS);
225: }

227: /*@
228:   PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials (the space is spanned
229:   by polynomials whose degree in each variable is bounded by the given order), as opposed to polynomials (the space is
230:   spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).

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

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

239:   Level: intermediate

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

251: /*@
252:   PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor polynomials (the space is spanned
253:   by polynomials whose degree in each variable is bounded by the given order), as opposed to polynomials (the space is
254:   spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).

256:   Input Parameter:
257: . sp     - the function space object

259:   Output Parameter:
260: . tensor - `PETSC_TRUE` for a tensor polynomial space, `PETSC_FALSE` for a polynomial space

262:   Level: intermediate

264: .seealso: `PetscSpace`, `PetscSpacePolynomialSetTensor()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
265: @*/
266: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
267: {
268:   PetscFunctionBegin;
271:   PetscTryMethod(sp, "PetscSpacePolynomialGetTensor_C", (PetscSpace, PetscBool *), (sp, tensor));
272:   PetscFunctionReturn(PETSC_SUCCESS);
273: }

275: static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
276: {
277:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;

279:   PetscFunctionBegin;
280:   poly->tensor = tensor;
281:   PetscFunctionReturn(PETSC_SUCCESS);
282: }

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

288:   PetscFunctionBegin;
291:   *tensor = poly->tensor;
292:   PetscFunctionReturn(PETSC_SUCCESS);
293: }

295: static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp)
296: {
297:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
298:   PetscInt         Nc, dim, order;
299:   PetscBool        tensor;

301:   PetscFunctionBegin;
302:   PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
303:   PetscCall(PetscSpaceGetNumVariables(sp, &dim));
304:   PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
305:   PetscCall(PetscSpacePolynomialGetTensor(sp, &tensor));
306:   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);
307:   if (!poly->subspaces) PetscCall(PetscCalloc1(dim, &poly->subspaces));
308:   if (height <= dim) {
309:     if (!poly->subspaces[height - 1]) {
310:       PetscSpace  sub;
311:       const char *name;

313:       PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub));
314:       PetscCall(PetscObjectGetName((PetscObject)sp, &name));
315:       PetscCall(PetscObjectSetName((PetscObject)sub, name));
316:       PetscCall(PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL));
317:       PetscCall(PetscSpaceSetNumComponents(sub, Nc));
318:       PetscCall(PetscSpaceSetDegree(sub, order, PETSC_DETERMINE));
319:       PetscCall(PetscSpaceSetNumVariables(sub, dim - height));
320:       PetscCall(PetscSpacePolynomialSetTensor(sub, tensor));
321:       PetscCall(PetscSpaceSetUp(sub));
322:       poly->subspaces[height - 1] = sub;
323:     }
324:     *subsp = poly->subspaces[height - 1];
325:   } else {
326:     *subsp = NULL;
327:   }
328:   PetscFunctionReturn(PETSC_SUCCESS);
329: }

331: static PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
332: {
333:   PetscFunctionBegin;
334:   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_Polynomial;
335:   sp->ops->setup             = PetscSpaceSetUp_Polynomial;
336:   sp->ops->view              = PetscSpaceView_Polynomial;
337:   sp->ops->destroy           = PetscSpaceDestroy_Polynomial;
338:   sp->ops->getdimension      = PetscSpaceGetDimension_Polynomial;
339:   sp->ops->evaluate          = PetscSpaceEvaluate_Polynomial;
340:   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial;
341:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial));
342:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial));
343:   PetscFunctionReturn(PETSC_SUCCESS);
344: }

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

350:   Level: intermediate

352: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
353: M*/

355: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
356: {
357:   PetscSpace_Poly *poly;

359:   PetscFunctionBegin;
361:   PetscCall(PetscNew(&poly));
362:   sp->data = poly;

364:   poly->tensor    = PETSC_FALSE;
365:   poly->subspaces = NULL;

367:   PetscCall(PetscSpaceInitialize_Polynomial(sp));
368:   PetscFunctionReturn(PETSC_SUCCESS);
369: }