Actual source code: spaceptrimmed.c
1: #include <petsc/private/petscfeimpl.h>
3: static PetscErrorCode PetscSpaceSetFromOptions_Ptrimmed(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
4: {
5: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data;
7: PetscOptionsHead(PetscOptionsObject,"PetscSpace polynomial options");
8: PetscOptionsInt("-petscspace_ptrimmed_form_degree", "form degree of trimmed space", "PetscSpacePTrimmedSetFormDegree", pt->formDegree, &(pt->formDegree), NULL);
9: PetscOptionsTail();
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 %D%s-forms of degree %D (P-%D/\\%D)\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) {
45: PetscSpaceDestroy(&pt->subspaces[d]);
46: }
47: }
48: PetscFree(pt->subspaces);
49: PetscFree(pt);
50: return 0;
51: }
53: static PetscErrorCode PetscSpaceSetUp_Ptrimmed(PetscSpace sp)
54: {
55: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data;
56: PetscInt Nf;
58: if (pt->setupCalled) return 0;
60: PetscDTBinomialInt(sp->Nv, PetscAbsInt(pt->formDegree), &Nf);
61: if (sp->Nc == PETSC_DETERMINE) {
62: sp->Nc = Nf;
63: }
65: if (sp->Nc != Nf) {
66: PetscSpace subsp;
67: PetscInt nCopies = sp->Nc / Nf;
68: PetscInt Nv, deg, maxDeg;
69: PetscInt formDegree = pt->formDegree;
70: const char *prefix;
71: const char *name;
72: char subname[PETSC_MAX_PATH_LEN];
74: PetscSpaceSetType(sp, PETSCSPACESUM);
75: PetscSpaceSumSetConcatenate(sp, PETSC_TRUE);
76: PetscSpaceSumSetNumSubspaces(sp, nCopies);
77: PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp);
78: PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix);
79: PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix);
80: PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_");
81: if (((PetscObject)sp)->name) {
82: PetscObjectGetName((PetscObject)sp, &name);
83: PetscSNPrintf(subname, PETSC_MAX_PATH_LEN-1, "%s sum component", name);
84: PetscObjectSetName((PetscObject)subsp, subname);
85: } else {
86: PetscObjectSetName((PetscObject)subsp, "sum component");
87: }
88: PetscSpaceSetType(subsp, PETSCSPACEPTRIMMED);
89: PetscSpaceGetNumVariables(sp, &Nv);
90: PetscSpaceSetNumVariables(subsp, Nv);
91: PetscSpaceSetNumComponents(subsp, Nf);
92: PetscSpaceGetDegree(sp, °, &maxDeg);
93: PetscSpaceSetDegree(subsp, deg, maxDeg);
94: PetscSpacePTrimmedSetFormDegree(subsp, formDegree);
95: PetscSpaceSetUp(subsp);
96: for (PetscInt i = 0; i < nCopies; i++) {
97: PetscSpaceSumSetSubspace(sp, i, subsp);
98: }
99: PetscSpaceDestroy(&subsp);
100: PetscSpaceSetUp(sp);
101: return 0;
102: }
103: if (sp->degree == PETSC_DEFAULT) {
104: sp->degree = 0;
105: } else if (sp->degree < 0) {
106: SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Invalid negative degree %D", sp->degree);
107: }
108: sp->maxDegree = (pt->formDegree == 0 || PetscAbsInt(pt->formDegree) == sp->Nv) ? sp->degree : sp->degree + 1;
109: if (pt->formDegree == 0 || PetscAbsInt(pt->formDegree) == sp->Nv) {
110: // Convert to regular polynomial space
111: PetscSpaceSetType(sp, PETSCSPACEPOLYNOMIAL);
112: PetscSpaceSetUp(sp);
113: return 0;
114: }
115: pt->setupCalled = PETSC_TRUE;
116: return 0;
117: }
119: static PetscErrorCode PetscSpaceGetDimension_Ptrimmed(PetscSpace sp, PetscInt *dim)
120: {
121: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data;
122: PetscInt f;
123: PetscInt Nf;
125: f = pt->formDegree;
126: // For PetscSpace, degree refers to the largest complete polynomial degree contained in the space which
127: // is equal to the index of a P trimmed space only for 0-forms: otherwise, the index is degree + 1
128: PetscDTPTrimmedSize(sp->Nv, f == 0 ? sp->degree : sp->degree + 1, pt->formDegree, dim);
129: PetscDTBinomialInt(sp->Nv, PetscAbsInt(pt->formDegree), &Nf);
130: *dim *= (sp->Nc / Nf);
131: return 0;
132: }
134: /*
135: p in [0, npoints), i in [0, pdim), c in [0, Nc)
137: B[p][i][c] = B[p][i_scalar][c][c]
138: */
139: static PetscErrorCode PetscSpaceEvaluate_Ptrimmed(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
140: {
141: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data;
142: DM dm = sp->dm;
143: PetscInt jet, degree, Nf, Ncopies, Njet;
144: PetscInt Nc = sp->Nc;
145: PetscInt f;
146: PetscInt dim = sp->Nv;
147: PetscReal *eval;
148: PetscInt Nb;
150: if (!pt->setupCalled) {
151: PetscSpaceSetUp(sp);
152: PetscSpaceEvaluate(sp, npoints, points, B, D, H);
153: return 0;
154: }
155: if (H) {
156: jet = 2;
157: } else if (D) {
158: jet = 1;
159: } else {
160: jet = 0;
161: }
162: f = pt->formDegree;
163: degree = f == 0 ? sp->degree : sp->degree + 1;
164: PetscDTBinomialInt(dim, PetscAbsInt(f), &Nf);
165: Ncopies = Nc / Nf;
167: PetscDTBinomialInt(dim + jet, dim, &Njet);
168: PetscDTPTrimmedSize(dim, degree, f, &Nb);
169: DMGetWorkArray(dm, Nb * Nf * Njet * npoints, MPIU_REAL, &eval);
170: PetscDTPTrimmedEvalJet(dim, npoints, points, degree, f, jet, eval);
171: if (B) {
172: PetscInt p_strl = Nf*Nb;
173: PetscInt b_strl = Nf;
174: PetscInt v_strl = 1;
176: PetscInt b_strr = Nf*Njet*npoints;
177: PetscInt v_strr = Njet*npoints;
178: PetscInt p_strr = 1;
180: for (PetscInt v = 0; v < Nf; v++) {
181: for (PetscInt b = 0; b < Nb; b++) {
182: for (PetscInt p = 0; p < npoints; p++) {
183: B[p*p_strl + b*b_strl + v*v_strl] = eval[b*b_strr + v*v_strr + p*p_strr];
184: }
185: }
186: }
187: }
188: if (D) {
189: PetscInt p_strl = dim*Nf*Nb;
190: PetscInt b_strl = dim*Nf;
191: PetscInt v_strl = dim;
192: PetscInt d_strl = 1;
194: PetscInt b_strr = Nf*Njet*npoints;
195: PetscInt v_strr = Njet*npoints;
196: PetscInt d_strr = npoints;
197: PetscInt p_strr = 1;
199: for (PetscInt v = 0; v < Nf; v++) {
200: for (PetscInt d = 0; d < dim; d++) {
201: for (PetscInt b = 0; b < Nb; b++) {
202: for (PetscInt p = 0; p < npoints; p++) {
203: 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];
204: }
205: }
206: }
207: }
208: }
209: if (H) {
210: PetscInt p_strl = dim*dim*Nf*Nb;
211: PetscInt b_strl = dim*dim*Nf;
212: PetscInt v_strl = dim*dim;
213: PetscInt d1_strl = dim;
214: PetscInt d2_strl = 1;
216: PetscInt b_strr = Nf*Njet*npoints;
217: PetscInt v_strr = Njet*npoints;
218: PetscInt j_strr = npoints;
219: PetscInt p_strr = 1;
221: PetscInt *derivs;
222: PetscCalloc1(dim, &derivs);
223: for (PetscInt d1 = 0; d1 < dim; d1++) {
224: for (PetscInt d2 = 0; d2 < dim; d2++) {
225: PetscInt j;
226: derivs[d1]++;
227: derivs[d2]++;
228: PetscDTGradedOrderToIndex(dim, derivs, &j);
229: derivs[d1]--;
230: derivs[d2]--;
231: for (PetscInt v = 0; v < Nf; v++) {
232: for (PetscInt b = 0; b < Nb; b++) {
233: for (PetscInt p = 0; p < npoints; p++) {
234: 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];
235: }
236: }
237: }
238: }
239: }
240: PetscFree(derivs);
241: }
242: DMRestoreWorkArray(dm, Nb * Nf * Njet * npoints, MPIU_REAL, &eval);
243: return 0;
244: }
246: /*@
247: PetscSpacePTrimmedSetFormDegree - Set the form degree of the trimmed polynomials.
249: Input Parameters:
250: + sp - the function space object
251: - formDegree - the form degree
253: Options Database:
254: . -petscspace_ptrimmed_form_degree <int> - The trimmed polynomial form degree
256: Level: intermediate
258: .seealso: PetscDTAltV, PetscDTPTrimmedEvalJet(), PetscSpacePTrimmedGetFormDegree()
259: @*/
260: PetscErrorCode PetscSpacePTrimmedSetFormDegree(PetscSpace sp, PetscInt formDegree)
261: {
263: PetscTryMethod(sp,"PetscSpacePTrimmedSetFormDegree_C",(PetscSpace,PetscInt),(sp,formDegree));
264: return 0;
265: }
267: /*@
268: PetscSpacePTrimmedGetFormDegree - Get the form degree of the trimmed polynomials.
270: Input Parameters:
271: . sp - the function space object
273: Output Parameters:
274: . formDegee - the form degree
276: Level: intermediate
278: .seealso: PetscDTAltV, PetscDTPTrimmedEvalJet(), PetscSpacePTrimmedSetFormDegree()
279: @*/
280: PetscErrorCode PetscSpacePTrimmedGetFormDegree(PetscSpace sp, PetscInt *formDegree)
281: {
284: PetscTryMethod(sp,"PetscSpacePTrimmedGetFormDegree_C",(PetscSpace,PetscInt*),(sp,formDegree));
285: return 0;
286: }
288: static PetscErrorCode PetscSpacePTrimmedSetFormDegree_Ptrimmed(PetscSpace sp, PetscInt formDegree)
289: {
290: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data;
292: pt->formDegree = formDegree;
293: return 0;
294: }
296: static PetscErrorCode PetscSpacePTrimmedGetFormDegree_Ptrimmed(PetscSpace sp, PetscInt *formDegree)
297: {
298: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data;
302: *formDegree = pt->formDegree;
303: return 0;
304: }
306: static PetscErrorCode PetscSpaceGetHeightSubspace_Ptrimmed(PetscSpace sp, PetscInt height, PetscSpace *subsp)
307: {
308: PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data;
309: PetscInt dim;
311: PetscSpaceGetNumVariables(sp, &dim);
313: if (!pt->subspaces) PetscCalloc1(dim, &(pt->subspaces));
314: if ((dim - height) <= PetscAbsInt(pt->formDegree)) {
315: if (!pt->subspaces[height-1]) {
316: PetscInt Nc, degree, Nf, Ncopies, Nfsub;
317: PetscSpace sub;
318: const char *name;
320: PetscSpaceGetNumComponents(sp, &Nc);
321: PetscDTBinomialInt(dim, PetscAbsInt(pt->formDegree), &Nf);
322: PetscDTBinomialInt((dim-height), PetscAbsInt(pt->formDegree), &Nfsub);
323: Ncopies = Nf / Nc;
324: PetscSpaceGetDegree(sp, °ree, NULL);
326: PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
327: PetscObjectGetName((PetscObject) sp, &name);
328: PetscObjectSetName((PetscObject) sub, name);
329: PetscSpaceSetType(sub, PETSCSPACEPTRIMMED);
330: PetscSpaceSetNumComponents(sub, Nfsub * Ncopies);
331: PetscSpaceSetDegree(sub, degree, PETSC_DETERMINE);
332: PetscSpaceSetNumVariables(sub, dim-height);
333: PetscSpacePTrimmedSetFormDegree(sub, pt->formDegree);
334: PetscSpaceSetUp(sub);
335: pt->subspaces[height-1] = sub;
336: }
337: *subsp = pt->subspaces[height-1];
338: } else {
339: *subsp = NULL;
340: }
341: return 0;
342: }
344: static PetscErrorCode PetscSpaceInitialize_Ptrimmed(PetscSpace sp)
345: {
346: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePTrimmedGetFormDegree_C", PetscSpacePTrimmedGetFormDegree_Ptrimmed);
347: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePTrimmedSetFormDegree_C", PetscSpacePTrimmedSetFormDegree_Ptrimmed);
348: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Ptrimmed;
349: sp->ops->setup = PetscSpaceSetUp_Ptrimmed;
350: sp->ops->view = PetscSpaceView_Ptrimmed;
351: sp->ops->destroy = PetscSpaceDestroy_Ptrimmed;
352: sp->ops->getdimension = PetscSpaceGetDimension_Ptrimmed;
353: sp->ops->evaluate = PetscSpaceEvaluate_Ptrimmed;
354: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Ptrimmed;
355: return 0;
356: }
358: /*MC
359: PETSCSPACEPTRIMMED = "ptrimmed" - A PetscSpace object that encapsulates a trimmed polynomial space.
361: Level: intermediate
363: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType(), PetscDTPTrimmedEvalJet()
364: M*/
366: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Ptrimmed(PetscSpace sp)
367: {
368: PetscSpace_Ptrimmed *pt;
371: PetscNewLog(sp,&pt);
372: sp->data = pt;
374: pt->subspaces = NULL;
375: sp->Nc = PETSC_DETERMINE;
377: PetscSpaceInitialize_Ptrimmed(sp);
378: return 0;
379: }