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, °, &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, °ree, 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: }