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:   PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace polynomial options");
  8:   PetscOptionsInt("-petscspace_ptrimmed_form_degree", "form degree of trimmed space", "PetscSpacePTrimmedSetFormDegree", pt->formDegree, &(pt->formDegree), NULL);
  9:   PetscOptionsHeadEnd();
 10:   return 0;
 11: }

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

 18:   f       = pt->formDegree;
 19:   tdegree = f == 0 ? sp->degree : sp->degree + 1;
 20:   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));
 21:   return 0;
 22: }

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

 30:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
 31:   if (iascii) PetscSpacePTrimmedView_Ascii(sp, viewer);
 32:   return 0;
 33: }

 35: static PetscErrorCode PetscSpaceDestroy_Ptrimmed(PetscSpace sp)
 36: {
 37:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;

 39:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePTrimmedGetFormDegree_C", NULL);
 40:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePTrimmedSetFormDegree_C", NULL);
 41:   if (pt->subspaces) {
 42:     PetscInt d;

 44:     for (d = 0; d < sp->Nv; ++d) PetscSpaceDestroy(&pt->subspaces[d]);
 45:   }
 46:   PetscFree(pt->subspaces);
 47:   PetscFree(pt);
 48:   return 0;
 49: }

 51: static PetscErrorCode PetscSpaceSetUp_Ptrimmed(PetscSpace sp)
 52: {
 53:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
 54:   PetscInt             Nf;

 56:   if (pt->setupCalled) return 0;
 58:   PetscDTBinomialInt(sp->Nv, PetscAbsInt(pt->formDegree), &Nf);
 59:   if (sp->Nc == PETSC_DETERMINE) sp->Nc = Nf;
 61:   if (sp->Nc != Nf) {
 62:     PetscSpace  subsp;
 63:     PetscInt    nCopies = sp->Nc / Nf;
 64:     PetscInt    Nv, deg, maxDeg;
 65:     PetscInt    formDegree = pt->formDegree;
 66:     const char *prefix;
 67:     const char *name;
 68:     char        subname[PETSC_MAX_PATH_LEN];

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

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

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

126: /*
127:   p in [0, npoints), i in [0, pdim), c in [0, Nc)

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

142:   if (!pt->setupCalled) {
143:     PetscSpaceSetUp(sp);
144:     PetscSpaceEvaluate(sp, npoints, points, B, D, H);
145:     return 0;
146:   }
147:   if (H) {
148:     jet = 2;
149:   } else if (D) {
150:     jet = 1;
151:   } else {
152:     jet = 0;
153:   }
154:   f      = pt->formDegree;
155:   degree = f == 0 ? sp->degree : sp->degree + 1;
156:   PetscDTBinomialInt(dim, PetscAbsInt(f), &Nf);
157:   Ncopies = Nc / Nf;
159:   PetscDTBinomialInt(dim + jet, dim, &Njet);
160:   PetscDTPTrimmedSize(dim, degree, f, &Nb);
161:   DMGetWorkArray(dm, Nb * Nf * Njet * npoints, MPIU_REAL, &eval);
162:   PetscDTPTrimmedEvalJet(dim, npoints, points, degree, f, jet, eval);
163:   if (B) {
164:     PetscInt p_strl = Nf * Nb;
165:     PetscInt b_strl = Nf;
166:     PetscInt v_strl = 1;

168:     PetscInt b_strr = Nf * Njet * npoints;
169:     PetscInt v_strr = Njet * npoints;
170:     PetscInt p_strr = 1;

172:     for (PetscInt v = 0; v < Nf; v++) {
173:       for (PetscInt b = 0; b < Nb; b++) {
174:         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];
175:       }
176:     }
177:   }
178:   if (D) {
179:     PetscInt p_strl = dim * Nf * Nb;
180:     PetscInt b_strl = dim * Nf;
181:     PetscInt v_strl = dim;
182:     PetscInt d_strl = 1;

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

189:     for (PetscInt v = 0; v < Nf; v++) {
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 + v * v_strl + d * d_strl] = eval[b * b_strr + v * v_strr + (1 + d) * d_strr + p * p_strr];
193:         }
194:       }
195:     }
196:   }
197:   if (H) {
198:     PetscInt p_strl  = dim * dim * Nf * Nb;
199:     PetscInt b_strl  = dim * dim * Nf;
200:     PetscInt v_strl  = dim * dim;
201:     PetscInt d1_strl = dim;
202:     PetscInt d2_strl = 1;

204:     PetscInt b_strr = Nf * Njet * npoints;
205:     PetscInt v_strr = Njet * npoints;
206:     PetscInt j_strr = npoints;
207:     PetscInt p_strr = 1;

209:     PetscInt *derivs;
210:     PetscCalloc1(dim, &derivs);
211:     for (PetscInt d1 = 0; d1 < dim; d1++) {
212:       for (PetscInt d2 = 0; d2 < dim; d2++) {
213:         PetscInt j;
214:         derivs[d1]++;
215:         derivs[d2]++;
216:         PetscDTGradedOrderToIndex(dim, derivs, &j);
217:         derivs[d1]--;
218:         derivs[d2]--;
219:         for (PetscInt v = 0; v < Nf; v++) {
220:           for (PetscInt b = 0; b < Nb; b++) {
221:             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];
222:           }
223:         }
224:       }
225:     }
226:     PetscFree(derivs);
227:   }
228:   DMRestoreWorkArray(dm, Nb * Nf * Njet * npoints, MPIU_REAL, &eval);
229:   return 0;
230: }

232: /*@
233:   PetscSpacePTrimmedSetFormDegree - Set the form degree of the trimmed polynomials.

235:   Input Parameters:
236: + sp         - the function space object
237: - formDegree - the form degree

239:   Options Database Key:
240: . -petscspace_ptrimmed_form_degree <int> - The trimmed polynomial form degree

242:   Level: intermediate

244: .seealso: `PetscSpace`, `PetscDTAltV`, `PetscDTPTrimmedEvalJet()`, `PetscSpacePTrimmedGetFormDegree()`
245: @*/
246: PetscErrorCode PetscSpacePTrimmedSetFormDegree(PetscSpace sp, PetscInt formDegree)
247: {
249:   PetscTryMethod(sp, "PetscSpacePTrimmedSetFormDegree_C", (PetscSpace, PetscInt), (sp, formDegree));
250:   return 0;
251: }

253: /*@
254:   PetscSpacePTrimmedGetFormDegree - Get the form degree of the trimmed polynomials.

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

259:   Output Parameters:
260: . formDegee - the form degree

262:   Level: intermediate

264: .seealso: `PetscSpace`, `PetscDTAltV`, `PetscDTPTrimmedEvalJet()`, `PetscSpacePTrimmedSetFormDegree()`
265: @*/
266: PetscErrorCode PetscSpacePTrimmedGetFormDegree(PetscSpace sp, PetscInt *formDegree)
267: {
270:   PetscTryMethod(sp, "PetscSpacePTrimmedGetFormDegree_C", (PetscSpace, PetscInt *), (sp, formDegree));
271:   return 0;
272: }

274: static PetscErrorCode PetscSpacePTrimmedSetFormDegree_Ptrimmed(PetscSpace sp, PetscInt formDegree)
275: {
276:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;

278:   pt->formDegree = formDegree;
279:   return 0;
280: }

282: static PetscErrorCode PetscSpacePTrimmedGetFormDegree_Ptrimmed(PetscSpace sp, PetscInt *formDegree)
283: {
284:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;

288:   *formDegree = pt->formDegree;
289:   return 0;
290: }

292: static PetscErrorCode PetscSpaceGetHeightSubspace_Ptrimmed(PetscSpace sp, PetscInt height, PetscSpace *subsp)
293: {
294:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
295:   PetscInt             dim;

297:   PetscSpaceGetNumVariables(sp, &dim);
299:   if (!pt->subspaces) PetscCalloc1(dim, &(pt->subspaces));
300:   if ((dim - height) <= PetscAbsInt(pt->formDegree)) {
301:     if (!pt->subspaces[height - 1]) {
302:       PetscInt    Nc, degree, Nf, Ncopies, Nfsub;
303:       PetscSpace  sub;
304:       const char *name;

306:       PetscSpaceGetNumComponents(sp, &Nc);
307:       PetscDTBinomialInt(dim, PetscAbsInt(pt->formDegree), &Nf);
308:       PetscDTBinomialInt((dim - height), PetscAbsInt(pt->formDegree), &Nfsub);
309:       Ncopies = Nf / Nc;
310:       PetscSpaceGetDegree(sp, &degree, NULL);

312:       PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub);
313:       PetscObjectGetName((PetscObject)sp, &name);
314:       PetscObjectSetName((PetscObject)sub, name);
315:       PetscSpaceSetType(sub, PETSCSPACEPTRIMMED);
316:       PetscSpaceSetNumComponents(sub, Nfsub * Ncopies);
317:       PetscSpaceSetDegree(sub, degree, PETSC_DETERMINE);
318:       PetscSpaceSetNumVariables(sub, dim - height);
319:       PetscSpacePTrimmedSetFormDegree(sub, pt->formDegree);
320:       PetscSpaceSetUp(sub);
321:       pt->subspaces[height - 1] = sub;
322:     }
323:     *subsp = pt->subspaces[height - 1];
324:   } else {
325:     *subsp = NULL;
326:   }
327:   return 0;
328: }

330: static PetscErrorCode PetscSpaceInitialize_Ptrimmed(PetscSpace sp)
331: {
332:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePTrimmedGetFormDegree_C", PetscSpacePTrimmedGetFormDegree_Ptrimmed);
333:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePTrimmedSetFormDegree_C", PetscSpacePTrimmedSetFormDegree_Ptrimmed);
334:   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_Ptrimmed;
335:   sp->ops->setup             = PetscSpaceSetUp_Ptrimmed;
336:   sp->ops->view              = PetscSpaceView_Ptrimmed;
337:   sp->ops->destroy           = PetscSpaceDestroy_Ptrimmed;
338:   sp->ops->getdimension      = PetscSpaceGetDimension_Ptrimmed;
339:   sp->ops->evaluate          = PetscSpaceEvaluate_Ptrimmed;
340:   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Ptrimmed;
341:   return 0;
342: }

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

347:   Level: intermediate

349: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscDTPTrimmedEvalJet()`
350: M*/

352: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Ptrimmed(PetscSpace sp)
353: {
354:   PetscSpace_Ptrimmed *pt;

357:   PetscNew(&pt);
358:   sp->data = pt;

360:   pt->subspaces = NULL;
361:   sp->Nc        = PETSC_DETERMINE;

363:   PetscSpaceInitialize_Ptrimmed(sp);
364:   return 0;
365: }