Actual source code: spacepoly.c

  1: #include <petsc/private/petscfeimpl.h>

  3: const char *const PetscSpacePolynomialTypes[] = {"P", "PMINUS_HDIV", "PMINUS_HCURL", "PetscSpacePolynomialType", "PETSCSPACE_POLYNOMIALTYPE_", NULL};

  5: static PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
  6: {
  7:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
  8:   PetscErrorCode   ierr;

 11:   PetscOptionsHead(PetscOptionsObject,"PetscSpace polynomial options");
 12:   PetscOptionsBool("-petscspace_poly_sym", "Use only symmetric polynomials", "PetscSpacePolynomialSetSymmetric", poly->symmetric, &poly->symmetric, NULL);
 13:   PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);
 14:   PetscOptionsEnum("-petscspace_poly_type", "Type of polynomial space", "PetscSpacePolynomialSetType", PetscSpacePolynomialTypes, (PetscEnum)poly->ptype, (PetscEnum*)&poly->ptype, NULL);
 15:   PetscOptionsTail();
 16:   return(0);
 17: }

 19: static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v)
 20: {
 21:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
 22:   PetscErrorCode   ierr;

 25:   PetscViewerASCIIPrintf(v, "%s%s%s space of degree %D\n", poly->ptype ? PetscSpacePolynomialTypes[poly->ptype] : "", poly->ptype ? " " : "", poly->tensor ? "Tensor polynomial" : "Polynomial", sp->degree);
 26:   return(0);
 27: }

 29: static PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
 30: {
 31:   PetscBool      iascii;

 37:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
 38:   if (iascii) {PetscSpacePolynomialView_Ascii(sp, viewer);}
 39:   return(0);
 40: }

 42: static PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
 43: {
 44:   PetscSpace_Poly *poly    = (PetscSpace_Poly *) sp->data;
 45:   PetscInt         ndegree = sp->degree+1;
 46:   PetscInt         deg;
 47:   PetscErrorCode   ierr;

 50:   if (poly->setupCalled) return(0);
 51:   PetscMalloc1(ndegree, &poly->degrees);
 52:   for (deg = 0; deg < ndegree; ++deg) poly->degrees[deg] = deg;
 53:   if (poly->tensor) {
 54:     sp->maxDegree = sp->degree + PetscMax(sp->Nv - 1,0);
 55:   } else {
 56:     sp->maxDegree = sp->degree;
 57:   }
 58:   poly->setupCalled = PETSC_TRUE;
 59:   return(0);
 60: }

 62: static PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
 63: {
 64:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
 65:   PetscErrorCode   ierr;

 68:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", NULL);
 69:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", NULL);
 70:   PetscFree(poly->degrees);
 71:   if (poly->subspaces) {
 72:     PetscInt d;

 74:     for (d = 0; d < sp->Nv; ++d) {
 75:       PetscSpaceDestroy(&poly->subspaces[d]);
 76:     }
 77:   }
 78:   PetscFree(poly->subspaces);
 79:   PetscFree(poly);
 80:   return(0);
 81: }

 83: /* We treat the space as a tensor product of scalar polynomial spaces, so the dimension is multiplied by Nc */
 84: static PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
 85: {
 86:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
 87:   PetscInt         deg  = sp->degree;
 88:   PetscInt         n    = sp->Nv, N, i;
 89:   PetscReal        D    = 1.0;

 92:   if ((poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HDIV) || (poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HCURL)) --deg;
 93:   if (poly->tensor) {
 94:     N = 1;
 95:     for (i = 0; i < n; ++i) N *= (deg+1);
 96:   } else {
 97:     for (i = 1; i <= n; ++i) {
 98:       D *= ((PetscReal) (deg+i))/i;
 99:     }
100:     N = (PetscInt) (D + 0.5);
101:   }
102:   if ((poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HDIV) || (poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HCURL)) {
103:     N *= sp->Nc + 1;
104:   } else {
105:     N *= sp->Nc;
106:   }
107:   *dim = N;
108:   return(0);
109: }

111: /*
112:   LatticePoint_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to 'sum'.

114:   Input Parameters:
115: + len - The length of the tuple
116: . sum - The sum of all entries in the tuple
117: - ind - The current multi-index of the tuple, initialized to the 0 tuple

119:   Output Parameter:
120: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
121: . tup - A tuple of len integers addig to sum

123:   Level: developer

125: .seealso:
126: */
127: static PetscErrorCode LatticePoint_Internal(PetscInt len, PetscInt sum, PetscInt ind[], PetscInt tup[])
128: {
129:   PetscInt       i;

133:   if (len == 1) {
134:     ind[0] = -1;
135:     tup[0] = sum;
136:   } else if (sum == 0) {
137:     for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
138:   } else {
139:     tup[0] = sum - ind[0];
140:     LatticePoint_Internal(len-1, ind[0], &ind[1], &tup[1]);
141:     if (ind[1] < 0) {
142:       if (ind[0] == sum) {ind[0] = -1;}
143:       else               {ind[1] = 0; ++ind[0];}
144:     }
145:   }
146:   return(0);
147: }

149: /*
150:   TensorPoint_Internal - Returns all tuples of size 'len' with nonnegative integers that are less than 'max'.

152:   Input Parameters:
153: + len - The length of the tuple
154: . max - The max for all entries in the tuple
155: - ind - The current multi-index of the tuple, initialized to the 0 tuple

157:   Output Parameter:
158: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
159: . tup - A tuple of len integers less than max

161:   Level: developer

163: .seealso:
164: */
165: static PetscErrorCode TensorPoint_Internal(PetscInt len, PetscInt max, PetscInt ind[], PetscInt tup[])
166: {
167:   PetscInt       i;

171:   if (len == 1) {
172:     tup[0] = ind[0]++;
173:     ind[0] = ind[0] >= max ? -1 : ind[0];
174:   } else if (max == 0) {
175:     for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
176:   } else {
177:     tup[0] = ind[0];
178:     TensorPoint_Internal(len-1, max, &ind[1], &tup[1]);
179:     if (ind[1] < 0) {
180:       ind[1] = 0;
181:       if (ind[0] == max-1) {ind[0] = -1;}
182:       else                 {++ind[0];}
183:     }
184:   }
185:   return(0);
186: }

188: /*
189:   p in [0, npoints), i in [0, pdim), c in [0, Nc)

191:   B[p][i][c] = B[p][i_scalar][c][c]
192: */
193: static PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
194: {
195:   const PetscInt eps[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}};
196:   PetscSpace_Poly *poly    = (PetscSpace_Poly *) sp->data;
197:   DM               dm      = sp->dm;
198:   PetscInt         Nc      = sp->Nc;
199:   PetscInt         ndegree = sp->degree+1;
200:   PetscInt        *degrees = poly->degrees;
201:   PetscInt         dim     = sp->Nv;
202:   PetscReal       *lpoints, *tmp, *LB, *LD, *LH;
203:   PetscInt        *ind, *tup;
204:   PetscInt         c, pdim, pdimRed, d, e, der, der2, i, p, deg, o;
205:   PetscErrorCode   ierr;

208:   PetscSpaceGetDimension(sp, &pdim);
209:   pdim /= Nc;
210:   DMGetWorkArray(dm, npoints, MPIU_REAL, &lpoints);
211:   DMGetWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);
212:   if (B || D || H) {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);}
213:   if (D || H)      {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);}
214:   if (H)           {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);}
215:   for (d = 0; d < dim; ++d) {
216:     for (p = 0; p < npoints; ++p) {
217:       lpoints[p] = points[p*dim+d];
218:     }
219:     PetscDTLegendreEval(npoints, lpoints, ndegree, degrees, tmp, &tmp[1*npoints*ndegree], &tmp[2*npoints*ndegree]);
220:     /* LB, LD, LH (ndegree * dim x npoints) */
221:     for (deg = 0; deg < ndegree; ++deg) {
222:       for (p = 0; p < npoints; ++p) {
223:         if (B || D || H) LB[(deg*dim + d)*npoints + p] = tmp[(0*npoints + p)*ndegree+deg];
224:         if (D || H)      LD[(deg*dim + d)*npoints + p] = tmp[(1*npoints + p)*ndegree+deg];
225:         if (H)           LH[(deg*dim + d)*npoints + p] = tmp[(2*npoints + p)*ndegree+deg];
226:       }
227:     }
228:   }
229:   /* Multiply by A (pdim x ndegree * dim) */
230:   PetscMalloc2(dim,&ind,dim,&tup);
231:   if (B) {
232:     PetscInt topDegree = sp->degree;

234:     /* B (npoints x pdim x Nc) */
235:     PetscArrayzero(B, npoints*pdim*Nc*Nc);
236:     if ((poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HDIV) || (poly->ptype == PETSCSPACE_POLYNOMIALTYPE_PMINUS_HCURL)) topDegree--;
237:     /* Make complete space portion */
238:     if (poly->tensor) {
239:       if (poly->ptype != PETSCSPACE_POLYNOMIALTYPE_P) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Tensor spaces not supported for P^- spaces (%s)", PetscSpacePolynomialTypes[poly->ptype]);
240:       i = 0;
241:       PetscArrayzero(ind, dim);
242:       while (ind[0] >= 0) {
243:         TensorPoint_Internal(dim, sp->degree+1, ind, tup);
244:         for (p = 0; p < npoints; ++p) {
245:           B[(p*pdim + i)*Nc*Nc] = 1.0;
246:           for (d = 0; d < dim; ++d) {
247:             B[(p*pdim + i)*Nc*Nc] *= LB[(tup[d]*dim + d)*npoints + p];
248:           }
249:         }
250:         ++i;
251:       }
252:     } else {
253:       i = 0;
254:       for (o = 0; o <= topDegree; ++o) {
255:         PetscArrayzero(ind, dim);
256:         while (ind[0] >= 0) {
257:           LatticePoint_Internal(dim, o, ind, tup);
258:           for (p = 0; p < npoints; ++p) {
259:             B[(p*pdim + i)*Nc*Nc] = 1.0;
260:             for (d = 0; d < dim; ++d) {
261:               B[(p*pdim + i)*Nc*Nc] *= LB[(tup[d]*dim + d)*npoints + p];
262:             }
263:           }
264:           ++i;
265:         }
266:       }
267:     }
268:     pdimRed = i;
269:     /* Make direct sum basis for multicomponent space */
270:     for (p = 0; p < npoints; ++p) {
271:       for (i = 0; i < pdimRed; ++i) {
272:         for (c = 1; c < Nc; ++c) {
273:           B[(p*pdim*Nc + i*Nc + c)*Nc + c] = B[(p*pdim + i)*Nc*Nc];
274:         }
275:       }
276:     }
277:     /* Make homogeneous part */
278:     if (topDegree < sp->degree) {
279:       if (poly->tensor) {
280:       } else {
281:         i = pdimRed;
282:         PetscArrayzero(ind, dim);
283:         while (ind[0] >= 0) {
284:           LatticePoint_Internal(dim, topDegree, ind, tup);
285:           for (p = 0; p < npoints; ++p) {
286:             for (c = 0; c < Nc; ++c) {
287:               B[(p*pdim*Nc + i*Nc + c)*Nc + c] = 1.0;
288:               for (d = 0; d < dim; ++d) {
289:                 B[(p*pdim*Nc + i*Nc + c)*Nc + c] *= LB[(tup[d]*dim + d)*npoints + p];
290:               }
291:               switch (poly->ptype) {
292:               case PETSCSPACE_POLYNOMIALTYPE_PMINUS_HDIV:
293:                 B[(p*pdim*Nc + i*Nc + c)*Nc + c] *= LB[(c*dim + d)*npoints + p];break;
294:               case PETSCSPACE_POLYNOMIALTYPE_PMINUS_HCURL:
295:               {
296:                 PetscReal sum = 0.0;
297:                 for (d = 0; d < dim; ++d) for (e = 0; e < dim; ++e) sum += eps[c][d][e]*LB[(d*dim + d)*npoints + p];
298:                 B[(p*pdim*Nc + i*Nc + c)*Nc + c] *= sum;
299:                 break;
300:               }
301:               default: SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Invalid polynomial type %s", PetscSpacePolynomialTypes[poly->ptype]);
302:               }
303:             }
304:           }
305:           ++i;
306:         }
307:       }
308:     }
309:   }
310:   if (D) {
311:     if (poly->ptype != PETSCSPACE_POLYNOMIALTYPE_P) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Derivatives not supported for P^- spaces (%s)", PetscSpacePolynomialTypes[poly->ptype]);
312:     /* D (npoints x pdim x Nc x dim) */
313:     PetscArrayzero(D, npoints*pdim*Nc*Nc*dim);
314:     if (poly->tensor) {
315:       i = 0;
316:       PetscArrayzero(ind, dim);
317:       while (ind[0] >= 0) {
318:         TensorPoint_Internal(dim, sp->degree+1, ind, tup);
319:         for (p = 0; p < npoints; ++p) {
320:           for (der = 0; der < dim; ++der) {
321:             D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0;
322:             for (d = 0; d < dim; ++d) {
323:               if (d == der) {
324:                 D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
325:               } else {
326:                 D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
327:               }
328:             }
329:           }
330:         }
331:         ++i;
332:       }
333:     } else {
334:       i = 0;
335:       for (o = 0; o <= sp->degree; ++o) {
336:         PetscArrayzero(ind, dim);
337:         while (ind[0] >= 0) {
338:           LatticePoint_Internal(dim, o, ind, tup);
339:           for (p = 0; p < npoints; ++p) {
340:             for (der = 0; der < dim; ++der) {
341:               D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0;
342:               for (d = 0; d < dim; ++d) {
343:                 if (d == der) {
344:                   D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
345:                 } else {
346:                   D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
347:                 }
348:               }
349:             }
350:           }
351:           ++i;
352:         }
353:       }
354:     }
355:     /* Make direct sum basis for multicomponent space */
356:     for (p = 0; p < npoints; ++p) {
357:       for (i = 0; i < pdim; ++i) {
358:         for (c = 1; c < Nc; ++c) {
359:           for (d = 0; d < dim; ++d) {
360:             D[((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d] = D[(p*pdim + i)*Nc*Nc*dim + d];
361:           }
362:         }
363:       }
364:     }
365:   }
366:   if (H) {
367:     if (poly->ptype != PETSCSPACE_POLYNOMIALTYPE_P) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Hessians not supported for P^- spaces (%s)", PetscSpacePolynomialTypes[poly->ptype]);
368:     /* H (npoints x pdim x Nc x Nc x dim x dim) */
369:     PetscArrayzero(H, npoints*pdim*Nc*Nc*dim*dim);
370:     if (poly->tensor) {
371:       i = 0;
372:       PetscArrayzero(ind, dim);
373:       while (ind[0] >= 0) {
374:         TensorPoint_Internal(dim, sp->degree+1, ind, tup);
375:         for (p = 0; p < npoints; ++p) {
376:           for (der = 0; der < dim; ++der) {
377:             H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] = 1.0;
378:             for (d = 0; d < dim; ++d) {
379:               if (d == der) {
380:                 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LH[(tup[d]*dim + d)*npoints + p];
381:               } else {
382:                 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
383:               }
384:             }
385:             for (der2 = der + 1; der2 < dim; ++der2) {
386:               H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0;
387:               for (d = 0; d < dim; ++d) {
388:                 if (d == der || d == der2) {
389:                   H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p];
390:                 } else {
391:                   H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p];
392:                 }
393:               }
394:               H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2];
395:             }
396:           }
397:         }
398:         ++i;
399:       }
400:     } else {
401:       i = 0;
402:       for (o = 0; o <= sp->degree; ++o) {
403:         PetscArrayzero(ind, dim);
404:         while (ind[0] >= 0) {
405:           LatticePoint_Internal(dim, o, ind, tup);
406:           for (p = 0; p < npoints; ++p) {
407:             for (der = 0; der < dim; ++der) {
408:               H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] = 1.0;
409:               for (d = 0; d < dim; ++d) {
410:                 if (d == der) {
411:                   H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LH[(tup[d]*dim + d)*npoints + p];
412:                 } else {
413:                   H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
414:                 }
415:               }
416:               for (der2 = der + 1; der2 < dim; ++der2) {
417:                 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0;
418:                 for (d = 0; d < dim; ++d) {
419:                   if (d == der || d == der2) {
420:                     H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p];
421:                   } else {
422:                     H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p];
423:                   }
424:                 }
425:                 H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2];
426:               }
427:             }
428:           }
429:           ++i;
430:         }
431:       }
432:     }
433:     /* Make direct sum basis for multicomponent space */
434:     for (p = 0; p < npoints; ++p) {
435:       for (i = 0; i < pdim; ++i) {
436:         for (c = 1; c < Nc; ++c) {
437:           for (d = 0; d < dim; ++d) {
438:             for (e = 0; e < dim; ++e) {
439:               H[(((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d)*dim + e] = H[((p*pdim + i)*Nc*Nc*dim + d)*dim + e];
440:             }
441:           }
442:         }
443:       }
444:     }
445:   }
446:   PetscFree2(ind,tup);
447:   if (H)           {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);}
448:   if (D || H)      {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);}
449:   if (B || D || H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);}
450:   DMRestoreWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);
451:   DMRestoreWorkArray(dm, npoints, MPIU_REAL, &lpoints);
452:   return(0);
453: }

455: /*@
456:   PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials (the space is spanned
457:   by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
458:   spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).

460:   Input Parameters:
461: + sp     - the function space object
462: - tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space

464:   Options Database:
465: . -petscspace_poly_tensor <bool> - Whether to use tensor product polynomials in higher dimension

467:   Level: intermediate

469: .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
470: @*/
471: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
472: {

477:   PetscTryMethod(sp,"PetscSpacePolynomialSetTensor_C",(PetscSpace,PetscBool),(sp,tensor));
478:   return(0);
479: }

481: /*@
482:   PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor polynomials (the space is spanned
483:   by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
484:   spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).

486:   Input Parameters:
487: . sp     - the function space object

489:   Output Parameters:
490: . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space

492:   Level: intermediate

494: .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
495: @*/
496: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
497: {

503:   PetscTryMethod(sp,"PetscSpacePolynomialGetTensor_C",(PetscSpace,PetscBool*),(sp,tensor));
504:   return(0);
505: }

507: static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
508: {
509:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

512:   poly->tensor = tensor;
513:   return(0);
514: }

516: static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor)
517: {
518:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

523:   *tensor = poly->tensor;
524:   return(0);
525: }

527: static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp)
528: {
529:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
530:   PetscInt         Nc, dim, order;
531:   PetscBool        tensor;
532:   PetscErrorCode   ierr;

535:   PetscSpaceGetNumComponents(sp, &Nc);
536:   PetscSpaceGetNumVariables(sp, &dim);
537:   PetscSpaceGetDegree(sp, &order, NULL);
538:   PetscSpacePolynomialGetTensor(sp, &tensor);
539:   if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim);
540:   if (!poly->subspaces) {PetscCalloc1(dim, &poly->subspaces);}
541:   if (height <= dim) {
542:     if (!poly->subspaces[height-1]) {
543:       PetscSpace  sub;
544:       const char *name;

546:       PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
547:       PetscObjectGetName((PetscObject) sp,  &name);
548:       PetscObjectSetName((PetscObject) sub,  name);
549:       PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL);
550:       PetscSpaceSetNumComponents(sub, Nc);
551:       PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
552:       PetscSpaceSetNumVariables(sub, dim-height);
553:       PetscSpacePolynomialSetTensor(sub, tensor);
554:       PetscSpaceSetUp(sub);
555:       poly->subspaces[height-1] = sub;
556:     }
557:     *subsp = poly->subspaces[height-1];
558:   } else {
559:     *subsp = NULL;
560:   }
561:   return(0);
562: }

564: static PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
565: {

569:   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_Polynomial;
570:   sp->ops->setup             = PetscSpaceSetUp_Polynomial;
571:   sp->ops->view              = PetscSpaceView_Polynomial;
572:   sp->ops->destroy           = PetscSpaceDestroy_Polynomial;
573:   sp->ops->getdimension      = PetscSpaceGetDimension_Polynomial;
574:   sp->ops->evaluate          = PetscSpaceEvaluate_Polynomial;
575:   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial;
576:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial);
577:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial);
578:   return(0);
579: }

581: /*MC
582:   PETSCSPACEPOLYNOMIAL = "poly" - A PetscSpace object that encapsulates a polynomial space, e.g. P1 is the space of
583:   linear polynomials. The space is replicated for each component.

585:   Level: intermediate

587: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
588: M*/

590: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
591: {
592:   PetscSpace_Poly *poly;
593:   PetscErrorCode   ierr;

597:   PetscNewLog(sp,&poly);
598:   sp->data = poly;

600:   poly->symmetric    = PETSC_FALSE;
601:   poly->tensor       = PETSC_FALSE;
602:   poly->degrees      = NULL;
603:   poly->subspaces    = NULL;

605:   PetscSpaceInitialize_Polynomial(sp);
606:   return(0);
607: }

609: /*@
610:   PetscSpacePolynomialSetSymmetric - Set whether a function space is a space of symmetric polynomials

612:   Input Parameters:
613: + sp  - the function space object
614: - sym - flag for symmetric polynomials

616:   Options Database:
617: . -petscspace_poly_sym <bool> - Whether to use symmetric polynomials

619:   Level: intermediate

621: .seealso: PetscSpacePolynomialGetSymmetric(), PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
622: @*/
623: PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym)
624: {
625:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

629:   poly->symmetric = sym;
630:   return(0);
631: }

633: /*@
634:   PetscSpacePolynomialGetSymmetric - Get whether a function space is a space of symmetric polynomials

636:   Input Parameter:
637: . sp  - the function space object

639:   Output Parameter:
640: . sym - flag for symmetric polynomials

642:   Level: intermediate

644: .seealso: PetscSpacePolynomialSetSymmetric(), PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
645: @*/
646: PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym)
647: {
648:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

653:   *sym = poly->symmetric;
654:   return(0);
655: }