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 isascii;
30: PetscFunctionBegin;
33: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
34: if (isascii) 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: for (PetscInt d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&pt->subspaces[d]));
47: }
48: PetscCall(PetscFree(pt->subspaces));
49: PetscCall(PetscFree(pt));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode PetscSpaceSetUp_Ptrimmed(PetscSpace sp)
54: {
55: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
56: PetscInt Nf;
58: PetscFunctionBegin;
59: if (pt->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
60: 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);
61: PetscCall(PetscDTBinomialInt(sp->Nv, PetscAbsInt(pt->formDegree), &Nf));
62: if (sp->Nc == PETSC_DETERMINE) sp->Nc = Nf;
63: 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);
64: if (sp->Nc != Nf) {
65: PetscSpace subsp;
66: PetscInt nCopies = sp->Nc / Nf;
67: PetscInt Nv, deg, maxDeg;
68: PetscInt formDegree = pt->formDegree;
69: const char *prefix;
70: const char *name;
71: char subname[PETSC_MAX_PATH_LEN];
73: PetscCall(PetscSpaceSetType(sp, PETSCSPACESUM));
74: PetscCall(PetscSpaceSumSetConcatenate(sp, PETSC_TRUE));
75: PetscCall(PetscSpaceSumSetNumSubspaces(sp, nCopies));
76: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp));
77: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix));
78: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix));
79: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_"));
80: if (((PetscObject)sp)->name) {
81: PetscCall(PetscObjectGetName((PetscObject)sp, &name));
82: PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s sum component", name));
83: PetscCall(PetscObjectSetName((PetscObject)subsp, subname));
84: } else PetscCall(PetscObjectSetName((PetscObject)subsp, "sum component"));
85: PetscCall(PetscSpaceSetType(subsp, PETSCSPACEPTRIMMED));
86: PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
87: PetscCall(PetscSpaceSetNumVariables(subsp, Nv));
88: PetscCall(PetscSpaceSetNumComponents(subsp, Nf));
89: PetscCall(PetscSpaceGetDegree(sp, °, &maxDeg));
90: PetscCall(PetscSpaceSetDegree(subsp, deg, maxDeg));
91: PetscCall(PetscSpacePTrimmedSetFormDegree(subsp, formDegree));
92: PetscCall(PetscSpaceSetUp(subsp));
93: for (PetscInt i = 0; i < nCopies; i++) PetscCall(PetscSpaceSumSetSubspace(sp, i, subsp));
94: PetscCall(PetscSpaceDestroy(&subsp));
95: PetscCall(PetscSpaceSetUp(sp));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
98: if (sp->degree == PETSC_DEFAULT) sp->degree = 0;
99: else PetscCheck(sp->degree >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Invalid negative degree %" PetscInt_FMT, sp->degree);
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: PetscCall(PetscSpaceSetType(sp, PETSCSPACEPOLYNOMIAL));
104: PetscCall(PetscSpaceSetUp(sp));
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
107: pt->setupcalled = PETSC_TRUE;
108: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
118: f = pt->formDegree;
119: // For PetscSpace, degree refers to the largest complete polynomial degree contained in the space which
120: // is equal to the index of a P trimmed space only for 0-forms: otherwise, the index is degree + 1
121: PetscCall(PetscDTPTrimmedSize(sp->Nv, f == 0 ? sp->degree : sp->degree + 1, pt->formDegree, dim));
122: PetscCall(PetscDTBinomialInt(sp->Nv, PetscAbsInt(pt->formDegree), &Nf));
123: *dim *= (sp->Nc / Nf);
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: /*
128: p in [0, npoints), i in [0, pdim), c in [0, Nc)
130: B[p][i][c] = B[p][i_scalar][c][c]
131: */
132: static PetscErrorCode PetscSpaceEvaluate_Ptrimmed(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
133: {
134: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
135: DM dm = sp->dm;
136: PetscInt jet, degree, Nf, Ncopies, Njet;
137: PetscInt Nc = sp->Nc;
138: PetscInt f;
139: PetscInt dim = sp->Nv;
140: PetscReal *eval;
141: PetscInt Nb;
143: PetscFunctionBegin;
144: if (!pt->setupcalled) {
145: PetscCall(PetscSpaceSetUp(sp));
146: PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
149: if (H) {
150: jet = 2;
151: } else if (D) {
152: jet = 1;
153: } else {
154: jet = 0;
155: }
156: f = pt->formDegree;
157: degree = f == 0 ? sp->degree : sp->degree + 1;
158: PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(f), &Nf));
159: Ncopies = Nc / Nf;
160: PetscCheck(Ncopies == 1, PetscObjectComm((PetscObject)sp), PETSC_ERR_PLIB, "Multicopy spaces should have been converted to PETSCSPACESUM");
161: PetscCall(PetscDTBinomialInt(dim + jet, dim, &Njet));
162: PetscCall(PetscDTPTrimmedSize(dim, degree, f, &Nb));
163: PetscCall(DMGetWorkArray(dm, Nb * Nf * Njet * npoints, MPIU_REAL, &eval));
164: PetscCall(PetscDTPTrimmedEvalJet(dim, npoints, points, degree, f, jet, eval));
165: if (B) {
166: PetscInt p_strl = Nf * Nb;
167: PetscInt b_strl = Nf;
168: PetscInt v_strl = 1;
170: PetscInt b_strr = Nf * Njet * npoints;
171: PetscInt v_strr = Njet * npoints;
172: PetscInt p_strr = 1;
174: for (PetscInt v = 0; v < Nf; v++) {
175: for (PetscInt b = 0; b < Nb; b++) {
176: 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];
177: }
178: }
179: }
180: if (D) {
181: PetscInt p_strl = dim * Nf * Nb;
182: PetscInt b_strl = dim * Nf;
183: PetscInt v_strl = dim;
184: PetscInt d_strl = 1;
186: PetscInt b_strr = Nf * Njet * npoints;
187: PetscInt v_strr = Njet * npoints;
188: PetscInt d_strr = npoints;
189: PetscInt p_strr = 1;
191: for (PetscInt v = 0; v < Nf; v++) {
192: for (PetscInt d = 0; d < dim; d++) {
193: for (PetscInt b = 0; b < Nb; b++) {
194: 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];
195: }
196: }
197: }
198: }
199: if (H) {
200: PetscInt p_strl = dim * dim * Nf * Nb;
201: PetscInt b_strl = dim * dim * Nf;
202: PetscInt v_strl = dim * dim;
203: PetscInt d1_strl = dim;
204: PetscInt d2_strl = 1;
206: PetscInt b_strr = Nf * Njet * npoints;
207: PetscInt v_strr = Njet * npoints;
208: PetscInt j_strr = npoints;
209: PetscInt p_strr = 1;
211: PetscInt *derivs;
212: PetscCall(PetscCalloc1(dim, &derivs));
213: for (PetscInt d1 = 0; d1 < dim; d1++) {
214: for (PetscInt d2 = 0; d2 < dim; d2++) {
215: PetscInt j;
216: derivs[d1]++;
217: derivs[d2]++;
218: PetscCall(PetscDTGradedOrderToIndex(dim, derivs, &j));
219: derivs[d1]--;
220: derivs[d2]--;
221: for (PetscInt v = 0; v < Nf; v++) {
222: for (PetscInt b = 0; b < Nb; b++) {
223: 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];
224: }
225: }
226: }
227: }
228: PetscCall(PetscFree(derivs));
229: }
230: PetscCall(DMRestoreWorkArray(dm, Nb * Nf * Njet * npoints, MPIU_REAL, &eval));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*@
235: PetscSpacePTrimmedSetFormDegree - Set the form degree of the trimmed polynomials.
237: Input Parameters:
238: + sp - the function space object
239: - formDegree - the form degree
241: Options Database Key:
242: . -petscspace_ptrimmed_form_degree degree - The trimmed polynomial form degree
244: Level: intermediate
246: .seealso: `PetscSpace`, `PetscDTAltV`, `PetscDTPTrimmedEvalJet()`, `PetscSpacePTrimmedGetFormDegree()`
247: @*/
248: PetscErrorCode PetscSpacePTrimmedSetFormDegree(PetscSpace sp, PetscInt formDegree)
249: {
250: PetscFunctionBegin;
252: PetscTryMethod(sp, "PetscSpacePTrimmedSetFormDegree_C", (PetscSpace, PetscInt), (sp, formDegree));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@
257: PetscSpacePTrimmedGetFormDegree - Get the form degree of the trimmed polynomials.
259: Input Parameter:
260: . sp - the function space object
262: Output Parameter:
263: . formDegree - the form degree
265: Level: intermediate
267: .seealso: `PetscSpace`, `PetscDTAltV`, `PetscDTPTrimmedEvalJet()`, `PetscSpacePTrimmedSetFormDegree()`
268: @*/
269: PetscErrorCode PetscSpacePTrimmedGetFormDegree(PetscSpace sp, PetscInt *formDegree)
270: {
271: PetscFunctionBegin;
273: PetscAssertPointer(formDegree, 2);
274: PetscTryMethod(sp, "PetscSpacePTrimmedGetFormDegree_C", (PetscSpace, PetscInt *), (sp, formDegree));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: static PetscErrorCode PetscSpacePTrimmedSetFormDegree_Ptrimmed(PetscSpace sp, PetscInt formDegree)
279: {
280: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
282: PetscFunctionBegin;
283: pt->formDegree = formDegree;
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: static PetscErrorCode PetscSpacePTrimmedGetFormDegree_Ptrimmed(PetscSpace sp, PetscInt *formDegree)
288: {
289: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
291: PetscFunctionBegin;
293: PetscAssertPointer(formDegree, 2);
294: *formDegree = pt->formDegree;
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: static PetscErrorCode PetscSpaceGetHeightSubspace_Ptrimmed(PetscSpace sp, PetscInt height, PetscSpace *subsp)
299: {
300: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *)sp->data;
301: PetscInt dim;
303: PetscFunctionBegin;
304: PetscCall(PetscSpaceGetNumVariables(sp, &dim));
305: 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);
306: if (!pt->subspaces) PetscCall(PetscCalloc1(dim, &pt->subspaces));
307: if ((dim - height) <= PetscAbsInt(pt->formDegree)) {
308: if (!pt->subspaces[height - 1]) {
309: PetscInt Nc, degree, Nf, Ncopies, Nfsub;
310: PetscSpace sub;
311: const char *name;
313: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
314: PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(pt->formDegree), &Nf));
315: PetscCall(PetscDTBinomialInt(dim - height, PetscAbsInt(pt->formDegree), &Nfsub));
316: Ncopies = Nf / Nc;
317: PetscCall(PetscSpaceGetDegree(sp, °ree, NULL));
319: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub));
320: PetscCall(PetscObjectGetName((PetscObject)sp, &name));
321: PetscCall(PetscObjectSetName((PetscObject)sub, name));
322: PetscCall(PetscSpaceSetType(sub, PETSCSPACEPTRIMMED));
323: PetscCall(PetscSpaceSetNumComponents(sub, Nfsub * Ncopies));
324: PetscCall(PetscSpaceSetDegree(sub, degree, PETSC_DETERMINE));
325: PetscCall(PetscSpaceSetNumVariables(sub, dim - height));
326: PetscCall(PetscSpacePTrimmedSetFormDegree(sub, pt->formDegree));
327: PetscCall(PetscSpaceSetUp(sub));
328: pt->subspaces[height - 1] = sub;
329: }
330: *subsp = pt->subspaces[height - 1];
331: } else {
332: *subsp = NULL;
333: }
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: static PetscErrorCode PetscSpaceInitialize_Ptrimmed(PetscSpace sp)
338: {
339: PetscFunctionBegin;
340: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePTrimmedGetFormDegree_C", PetscSpacePTrimmedGetFormDegree_Ptrimmed));
341: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePTrimmedSetFormDegree_C", PetscSpacePTrimmedSetFormDegree_Ptrimmed));
342: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Ptrimmed;
343: sp->ops->setup = PetscSpaceSetUp_Ptrimmed;
344: sp->ops->view = PetscSpaceView_Ptrimmed;
345: sp->ops->destroy = PetscSpaceDestroy_Ptrimmed;
346: sp->ops->getdimension = PetscSpaceGetDimension_Ptrimmed;
347: sp->ops->evaluate = PetscSpaceEvaluate_Ptrimmed;
348: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Ptrimmed;
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: /*MC
353: PETSCSPACEPTRIMMED = "ptrimmed" - A `PetscSpace` object that encapsulates a trimmed polynomial space.
355: Trimmed polynomial spaces are defined for $k$-forms, and are defined by
356: $
357: \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)],
358: $
359: 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.
361: Level: intermediate
363: Notes:
364: Trimmed polynomial spaces correspond to several common conformal approximation spaces in the de Rham complex:
366: In $H^1$ ($\sim k=0$), trimmed polynomial spaces are identical to the standard polynomial spaces, $\mathcal{P}_r^- \sim P_r$.
368: 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
369: $
370: \begin{cases}
371: [P_{r-1}(\mathbb{R}^2)]^2 \oplus \mathrm{rot}(\bf{x}) H_{r-1}(\mathbb{R}^2), & n = 2, \\
372: [P_{r-1}(\mathbb{R}^3)]^3 \oplus \bf{x} \times [H_{r-1}(\mathbb{R}^3)]^3, & n = 3.
373: \end{cases}
374: $
376: 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
377: $
378: [P_{r-1}(\mathbb{R}^n)]^n \oplus \bf{x} H_{r-1}(\mathbb{R}^n).
379: $
381: In $L_2$, ($\sim k=n$), trimmed polynomial spaces are identical to the standard polynomial spaces of one degree less, $\mathcal{P}_r^- \sim P_{r-1}$.
383: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscDTPTrimmedEvalJet()`
384: M*/
386: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Ptrimmed(PetscSpace sp)
387: {
388: PetscSpace_Ptrimmed *pt;
390: PetscFunctionBegin;
392: PetscCall(PetscNew(&pt));
393: sp->data = pt;
395: pt->subspaces = NULL;
396: sp->Nc = PETSC_DETERMINE;
398: PetscCall(PetscSpaceInitialize_Ptrimmed(sp));
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }