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, °, &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, °ree, 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: Level: intermediate
362: Developer Note:
363: Need a good easy to understand reference for trimmed poynomial spaces
365: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscDTPTrimmedEvalJet()`
366: M*/
368: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Ptrimmed(PetscSpace sp)
369: {
370: PetscSpace_Ptrimmed *pt;
372: PetscFunctionBegin;
374: PetscCall(PetscNew(&pt));
375: sp->data = pt;
377: pt->subspaces = NULL;
378: sp->Nc = PETSC_DETERMINE;
380: PetscCall(PetscSpaceInitialize_Ptrimmed(sp));
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }