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 iascii;

 30:   PetscFunctionBegin;
 33:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 34:   if (iascii) 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) {
101:     sp->degree = 0;
102:   } else if (sp->degree < 0) {
103:     SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Invalid negative degree %" PetscInt_FMT, sp->degree);
104:   }
105:   sp->maxDegree = (pt->formDegree == 0 || PetscAbsInt(pt->formDegree) == sp->Nv) ? sp->degree : sp->degree + 1;
106:   if (pt->formDegree == 0 || PetscAbsInt(pt->formDegree) == sp->Nv) {
107:     // Convert to regular polynomial space
108:     PetscCall(PetscSpaceSetType(sp, PETSCSPACEPOLYNOMIAL));
109:     PetscCall(PetscSpaceSetUp(sp));
110:     PetscFunctionReturn(PETSC_SUCCESS);
111:   }
112:   pt->setupCalled = PETSC_TRUE;
113:   PetscFunctionReturn(PETSC_SUCCESS);
114: }

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

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

132: /*
133:   p in [0, npoints), i in [0, pdim), c in [0, Nc)

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

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

175:     PetscInt b_strr = Nf * Njet * npoints;
176:     PetscInt v_strr = Njet * npoints;
177:     PetscInt p_strr = 1;

179:     for (PetscInt v = 0; v < Nf; v++) {
180:       for (PetscInt b = 0; b < Nb; b++) {
181:         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];
182:       }
183:     }
184:   }
185:   if (D) {
186:     PetscInt p_strl = dim * Nf * Nb;
187:     PetscInt b_strl = dim * Nf;
188:     PetscInt v_strl = dim;
189:     PetscInt d_strl = 1;

191:     PetscInt b_strr = Nf * Njet * npoints;
192:     PetscInt v_strr = Njet * npoints;
193:     PetscInt d_strr = npoints;
194:     PetscInt p_strr = 1;

196:     for (PetscInt v = 0; v < Nf; v++) {
197:       for (PetscInt d = 0; d < dim; d++) {
198:         for (PetscInt b = 0; b < Nb; b++) {
199:           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];
200:         }
201:       }
202:     }
203:   }
204:   if (H) {
205:     PetscInt p_strl  = dim * dim * Nf * Nb;
206:     PetscInt b_strl  = dim * dim * Nf;
207:     PetscInt v_strl  = dim * dim;
208:     PetscInt d1_strl = dim;
209:     PetscInt d2_strl = 1;

211:     PetscInt b_strr = Nf * Njet * npoints;
212:     PetscInt v_strr = Njet * npoints;
213:     PetscInt j_strr = npoints;
214:     PetscInt p_strr = 1;

216:     PetscInt *derivs;
217:     PetscCall(PetscCalloc1(dim, &derivs));
218:     for (PetscInt d1 = 0; d1 < dim; d1++) {
219:       for (PetscInt d2 = 0; d2 < dim; d2++) {
220:         PetscInt j;
221:         derivs[d1]++;
222:         derivs[d2]++;
223:         PetscCall(PetscDTGradedOrderToIndex(dim, derivs, &j));
224:         derivs[d1]--;
225:         derivs[d2]--;
226:         for (PetscInt v = 0; v < Nf; v++) {
227:           for (PetscInt b = 0; b < Nb; b++) {
228:             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];
229:           }
230:         }
231:       }
232:     }
233:     PetscCall(PetscFree(derivs));
234:   }
235:   PetscCall(DMRestoreWorkArray(dm, Nb * Nf * Njet * npoints, MPIU_REAL, &eval));
236:   PetscFunctionReturn(PETSC_SUCCESS);
237: }

239: /*@
240:   PetscSpacePTrimmedSetFormDegree - Set the form degree of the trimmed polynomials.

242:   Input Parameters:
243: + sp         - the function space object
244: - formDegree - the form degree

246:   Options Database Key:
247: . -petscspace_ptrimmed_form_degree <int> - The trimmed polynomial form degree

249:   Level: intermediate

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

261: /*@
262:   PetscSpacePTrimmedGetFormDegree - Get the form degree of the trimmed polynomials.

264:   Input Parameter:
265: . sp - the function space object

267:   Output Parameter:
268: . formDegree - the form degree

270:   Level: intermediate

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

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

287:   PetscFunctionBegin;
288:   pt->formDegree = formDegree;
289:   PetscFunctionReturn(PETSC_SUCCESS);
290: }

292: static PetscErrorCode PetscSpacePTrimmedGetFormDegree_Ptrimmed(PetscSpace sp, PetscInt *formDegree)
293: {
294:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;

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

303: static PetscErrorCode PetscSpaceGetHeightSubspace_Ptrimmed(PetscSpace sp, PetscInt height, PetscSpace *subsp)
304: {
305:   PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
306:   PetscInt             dim;

308:   PetscFunctionBegin;
309:   PetscCall(PetscSpaceGetNumVariables(sp, &dim));
310:   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);
311:   if (!pt->subspaces) PetscCall(PetscCalloc1(dim, &pt->subspaces));
312:   if ((dim - height) <= PetscAbsInt(pt->formDegree)) {
313:     if (!pt->subspaces[height - 1]) {
314:       PetscInt    Nc, degree, Nf, Ncopies, Nfsub;
315:       PetscSpace  sub;
316:       const char *name;

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

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

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

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

360:   Trimmed polynomial spaces are defined for $k$-forms, and are defined by
361:   $
362:   \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)],
363:   $
364:   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.

366:   Level: intermediate

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

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

373:   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
374:   $
375:     \begin{cases}
376:       [P_{r-1}(\mathbb{R}^2)]^2 \oplus \mathrm{rot}(\bf{x}) H_{r-1}(\mathbb{R}^2), & n = 2, \\
377:       [P_{r-1}(\mathbb{R}^3)]^3 \oplus \bf{x} \times [H_{r-1}(\mathbb{R}^3)]^3, & n = 3.
378:     \end{cases}
379:   $

381:   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
382:   $
383:     [P_{r-1}(\mathbb{R}^n)]^n \oplus \bf{x} H_{r-1}(\mathbb{R}^n).
384:   $

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

388: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscDTPTrimmedEvalJet()`
389: M*/

391: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Ptrimmed(PetscSpace sp)
392: {
393:   PetscSpace_Ptrimmed *pt;

395:   PetscFunctionBegin;
397:   PetscCall(PetscNew(&pt));
398:   sp->data = pt;

400:   pt->subspaces = NULL;
401:   sp->Nc        = PETSC_DETERMINE;

403:   PetscCall(PetscSpaceInitialize_Ptrimmed(sp));
404:   PetscFunctionReturn(PETSC_SUCCESS);
405: }