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:   PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace polynomial options");
  8:   PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);
  9:   PetscOptionsHeadEnd();
 10:   return 0;
 11: }

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

 17:   PetscViewerASCIIPrintf(v, "%s space of degree %" PetscInt_FMT "\n", poly->tensor ? "Tensor polynomial" : "Polynomial", sp->degree);
 18:   return 0;
 19: }

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

 27:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
 28:   if (iascii) PetscSpacePolynomialView_Ascii(sp, viewer);
 29:   return 0;
 30: }

 32: static PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
 33: {
 34:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;

 36:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", NULL);
 37:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialSetTensor_C", NULL);
 38:   if (poly->subspaces) {
 39:     PetscInt d;

 41:     for (d = 0; d < sp->Nv; ++d) PetscSpaceDestroy(&poly->subspaces[d]);
 42:   }
 43:   PetscFree(poly->subspaces);
 44:   PetscFree(poly);
 45:   return 0;
 46: }

 48: static PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
 49: {
 50:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;

 52:   if (poly->setupCalled) return 0;
 53:   if (sp->Nv <= 1) poly->tensor = PETSC_FALSE;
 54:   if (sp->Nc != 1) {
 55:     PetscInt    Nc     = sp->Nc;
 56:     PetscBool   tensor = poly->tensor;
 57:     PetscInt    Nv     = sp->Nv;
 58:     PetscInt    degree = sp->degree;
 59:     const char *prefix;
 60:     const char *name;
 61:     char        subname[PETSC_MAX_PATH_LEN];
 62:     PetscSpace  subsp;

 64:     PetscSpaceSetType(sp, PETSCSPACESUM);
 65:     PetscSpaceSumSetNumSubspaces(sp, Nc);
 66:     PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp);
 67:     PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix);
 68:     PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix);
 69:     PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_");
 70:     if (((PetscObject)sp)->name) {
 71:       PetscObjectGetName((PetscObject)sp, &name);
 72:       PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s sum component", name);
 73:       PetscObjectSetName((PetscObject)subsp, subname);
 74:     } else PetscObjectSetName((PetscObject)subsp, "sum component");
 75:     PetscSpaceSetType(subsp, PETSCSPACEPOLYNOMIAL);
 76:     PetscSpaceSetDegree(subsp, degree, PETSC_DETERMINE);
 77:     PetscSpaceSetNumComponents(subsp, 1);
 78:     PetscSpaceSetNumVariables(subsp, Nv);
 79:     PetscSpacePolynomialSetTensor(subsp, tensor);
 80:     PetscSpaceSetUp(subsp);
 81:     for (PetscInt i = 0; i < Nc; i++) PetscSpaceSumSetSubspace(sp, i, subsp);
 82:     PetscSpaceDestroy(&subsp);
 83:     PetscSpaceSetUp(sp);
 84:     return 0;
 85:   }
 86:   if (poly->tensor) {
 87:     sp->maxDegree = PETSC_DETERMINE;
 88:     PetscSpaceSetType(sp, PETSCSPACETENSOR);
 89:     PetscSpaceSetUp(sp);
 90:     return 0;
 91:   }
 93:   sp->maxDegree     = sp->degree;
 94:   poly->setupCalled = PETSC_TRUE;
 95:   return 0;
 96: }

 98: static PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
 99: {
100:   PetscInt deg = sp->degree;
101:   PetscInt n   = sp->Nv;

103:   PetscDTBinomialInt(n + deg, n, dim);
104:   *dim *= sp->Nc;
105:   return 0;
106: }

108: static PetscErrorCode CoordinateBasis(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt jet, PetscInt Njet, PetscReal pScalar[])
109: {
110:   PetscArrayzero(pScalar, (1 + dim) * Njet * npoints);
111:   for (PetscInt b = 0; b < 1 + dim; b++) {
112:     for (PetscInt j = 0; j < PetscMin(1 + dim, Njet); j++) {
113:       if (j == 0) {
114:         if (b == 0) {
115:           for (PetscInt pt = 0; pt < npoints; pt++) pScalar[b * Njet * npoints + j * npoints + pt] = 1.;
116:         } else {
117:           for (PetscInt pt = 0; pt < npoints; pt++) pScalar[b * Njet * npoints + j * npoints + pt] = points[pt * dim + (b - 1)];
118:         }
119:       } else if (j == b) {
120:         for (PetscInt pt = 0; pt < npoints; pt++) pScalar[b * Njet * npoints + j * npoints + pt] = 1.;
121:       }
122:     }
123:   }
124:   return 0;
125: }

127: static PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
128: {
129:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
130:   DM               dm   = sp->dm;
131:   PetscInt         dim  = sp->Nv;
132:   PetscInt         Nb, jet, Njet;
133:   PetscReal       *pScalar;

135:   if (!poly->setupCalled) {
136:     PetscSpaceSetUp(sp);
137:     PetscSpaceEvaluate(sp, npoints, points, B, D, H);
138:     return 0;
139:   }
141:   PetscDTBinomialInt(dim + sp->degree, dim, &Nb);
142:   if (H) {
143:     jet = 2;
144:   } else if (D) {
145:     jet = 1;
146:   } else {
147:     jet = 0;
148:   }
149:   PetscDTBinomialInt(dim + jet, dim, &Njet);
150:   DMGetWorkArray(dm, Nb * Njet * npoints, MPIU_REAL, &pScalar);
151:   // Why are we handling the case degree == 1 specially?  Because we don't want numerical noise when we evaluate hat
152:   // functions at the vertices of a simplex, which happens when we invert the Vandermonde matrix of the PKD basis.
153:   // We don't make any promise about which basis is used.
154:   if (sp->degree == 1) {
155:     CoordinateBasis(dim, npoints, points, jet, Njet, pScalar);
156:   } else {
157:     PetscDTPKDEvalJet(dim, npoints, points, sp->degree, jet, pScalar);
158:   }
159:   if (B) {
160:     PetscInt p_strl = Nb;
161:     PetscInt b_strl = 1;

163:     PetscInt b_strr = Njet * npoints;
164:     PetscInt p_strr = 1;

166:     PetscArrayzero(B, npoints * Nb);
167:     for (PetscInt b = 0; b < Nb; b++) {
168:       for (PetscInt p = 0; p < npoints; p++) B[p * p_strl + b * b_strl] = pScalar[b * b_strr + p * p_strr];
169:     }
170:   }
171:   if (D) {
172:     PetscInt p_strl = dim * Nb;
173:     PetscInt b_strl = dim;
174:     PetscInt d_strl = 1;

176:     PetscInt b_strr = Njet * npoints;
177:     PetscInt d_strr = npoints;
178:     PetscInt p_strr = 1;

180:     PetscArrayzero(D, npoints * Nb * dim);
181:     for (PetscInt d = 0; d < dim; d++) {
182:       for (PetscInt b = 0; b < Nb; b++) {
183:         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];
184:       }
185:     }
186:   }
187:   if (H) {
188:     PetscInt p_strl  = dim * dim * Nb;
189:     PetscInt b_strl  = dim * dim;
190:     PetscInt d1_strl = dim;
191:     PetscInt d2_strl = 1;

193:     PetscInt b_strr = Njet * npoints;
194:     PetscInt j_strr = npoints;
195:     PetscInt p_strr = 1;

197:     PetscInt *derivs;
198:     PetscCalloc1(dim, &derivs);
199:     PetscArrayzero(H, npoints * Nb * dim * dim);
200:     for (PetscInt d1 = 0; d1 < dim; d1++) {
201:       for (PetscInt d2 = 0; d2 < dim; d2++) {
202:         PetscInt j;
203:         derivs[d1]++;
204:         derivs[d2]++;
205:         PetscDTGradedOrderToIndex(dim, derivs, &j);
206:         derivs[d1]--;
207:         derivs[d2]--;
208:         for (PetscInt b = 0; b < Nb; b++) {
209:           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];
210:         }
211:       }
212:     }
213:     PetscFree(derivs);
214:   }
215:   DMRestoreWorkArray(dm, Nb * Njet * npoints, MPIU_REAL, &pScalar);
216:   return 0;
217: }

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

224:   Input Parameters:
225: + sp     - the function space object
226: - tensor - `PETSC_TRUE` for a tensor polynomial space, `PETSC_FALSE` for a polynomial space

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

231:   Level: intermediate

233: .seealso: `PetscSpace`, `PetscSpacePolynomialGetTensor()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
234: @*/
235: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
236: {
238:   PetscTryMethod(sp, "PetscSpacePolynomialSetTensor_C", (PetscSpace, PetscBool), (sp, tensor));
239:   return 0;
240: }

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

247:   Input Parameters:
248: . sp     - the function space object

250:   Output Parameters:
251: . tensor - `PETSC_TRUE` for a tensor polynomial space, `PETSC_FALSE` for a polynomial space

253:   Level: intermediate

255: .seealso: `PetscSpace`, `PetscSpacePolynomialSetTensor()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
256: @*/
257: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
258: {
261:   PetscTryMethod(sp, "PetscSpacePolynomialGetTensor_C", (PetscSpace, PetscBool *), (sp, tensor));
262:   return 0;
263: }

265: static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
266: {
267:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;

269:   poly->tensor = tensor;
270:   return 0;
271: }

273: static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor)
274: {
275:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;

279:   *tensor = poly->tensor;
280:   return 0;
281: }

283: static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp)
284: {
285:   PetscSpace_Poly *poly = (PetscSpace_Poly *)sp->data;
286:   PetscInt         Nc, dim, order;
287:   PetscBool        tensor;

289:   PetscSpaceGetNumComponents(sp, &Nc);
290:   PetscSpaceGetNumVariables(sp, &dim);
291:   PetscSpaceGetDegree(sp, &order, NULL);
292:   PetscSpacePolynomialGetTensor(sp, &tensor);
294:   if (!poly->subspaces) PetscCalloc1(dim, &poly->subspaces);
295:   if (height <= dim) {
296:     if (!poly->subspaces[height - 1]) {
297:       PetscSpace  sub;
298:       const char *name;

300:       PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub);
301:       PetscObjectGetName((PetscObject)sp, &name);
302:       PetscObjectSetName((PetscObject)sub, name);
303:       PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL);
304:       PetscSpaceSetNumComponents(sub, Nc);
305:       PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
306:       PetscSpaceSetNumVariables(sub, dim - height);
307:       PetscSpacePolynomialSetTensor(sub, tensor);
308:       PetscSpaceSetUp(sub);
309:       poly->subspaces[height - 1] = sub;
310:     }
311:     *subsp = poly->subspaces[height - 1];
312:   } else {
313:     *subsp = NULL;
314:   }
315:   return 0;
316: }

318: static PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
319: {
320:   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_Polynomial;
321:   sp->ops->setup             = PetscSpaceSetUp_Polynomial;
322:   sp->ops->view              = PetscSpaceView_Polynomial;
323:   sp->ops->destroy           = PetscSpaceDestroy_Polynomial;
324:   sp->ops->getdimension      = PetscSpaceGetDimension_Polynomial;
325:   sp->ops->evaluate          = PetscSpaceEvaluate_Polynomial;
326:   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial;
327:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial);
328:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial);
329:   return 0;
330: }

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

336:   Level: intermediate

338: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
339: M*/

341: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
342: {
343:   PetscSpace_Poly *poly;

346:   PetscNew(&poly);
347:   sp->data = poly;

349:   poly->tensor    = PETSC_FALSE;
350:   poly->subspaces = NULL;

352:   PetscSpaceInitialize_Polynomial(sp);
353:   return 0;
354: }