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, &deg, &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, &degree, 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: }