Actual source code: spacepoly.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1:  #include <petsc/private/petscfeimpl.h>

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

  9:   PetscOptionsHead(PetscOptionsObject,"PetscSpace polynomial options");
 10:   PetscOptionsBool("-petscspace_poly_sym", "Use only symmetric polynomials", "PetscSpacePolynomialSetSymmetric", poly->symmetric, &poly->symmetric, NULL);
 11:   PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);
 12:   PetscOptionsTail();
 13:   return(0);
 14: }

 16: static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v)
 17: {
 18:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
 19:   PetscErrorCode   ierr;

 22:   PetscViewerASCIIPrintf(v, "%s space of degree %D\n", poly->tensor ? "Tensor polynomial" : "Polynomial", sp->degree);
 23:   return(0);
 24: }

 26: PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
 27: {
 28:   PetscBool      iascii;

 34:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
 35:   if (iascii) {PetscSpacePolynomialView_Ascii(sp, viewer);}
 36:   return(0);
 37: }

 39: PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
 40: {
 41:   PetscSpace_Poly *poly    = (PetscSpace_Poly *) sp->data;
 42:   PetscInt         ndegree = sp->degree+1;
 43:   PetscInt         deg;
 44:   PetscErrorCode   ierr;

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

 59: PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
 60: {
 61:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
 62:   PetscErrorCode   ierr;

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

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

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

 89:   if (poly->tensor) {
 90:     *dim = 1;
 91:     for (i = 0; i < n; ++i) *dim *= (deg+1);
 92:   } else {
 93:     for (i = 1; i <= n; ++i) {
 94:       D *= ((PetscReal) (deg+i))/i;
 95:     }
 96:     *dim = (PetscInt) (D + 0.5);
 97:   }
 98:   *dim *= sp->Nc;
 99:   return(0);
100: }

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

105:   Input Parameters:
106: + len - The length of the tuple
107: . sum - The sum of all entries in the tuple
108: - ind - The current multi-index of the tuple, initialized to the 0 tuple

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

114:   Level: developer

116: .seealso:
117: */
118: static PetscErrorCode LatticePoint_Internal(PetscInt len, PetscInt sum, PetscInt ind[], PetscInt tup[])
119: {
120:   PetscInt       i;

124:   if (len == 1) {
125:     ind[0] = -1;
126:     tup[0] = sum;
127:   } else if (sum == 0) {
128:     for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
129:   } else {
130:     tup[0] = sum - ind[0];
131:     LatticePoint_Internal(len-1, ind[0], &ind[1], &tup[1]);
132:     if (ind[1] < 0) {
133:       if (ind[0] == sum) {ind[0] = -1;}
134:       else               {ind[1] = 0; ++ind[0];}
135:     }
136:   }
137:   return(0);
138: }

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

143:   Input Parameters:
144: + len - The length of the tuple
145: . max - The max for all entries in the tuple
146: - ind - The current multi-index of the tuple, initialized to the 0 tuple

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

152:   Level: developer

154: .seealso:
155: */
156: static PetscErrorCode TensorPoint_Internal(PetscInt len, PetscInt max, PetscInt ind[], PetscInt tup[])
157: {
158:   PetscInt       i;

162:   if (len == 1) {
163:     tup[0] = ind[0]++;
164:     ind[0] = ind[0] >= max ? -1 : ind[0];
165:   } else if (max == 0) {
166:     for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
167:   } else {
168:     tup[0] = ind[0];
169:     TensorPoint_Internal(len-1, max, &ind[1], &tup[1]);
170:     if (ind[1] < 0) {
171:       ind[1] = 0;
172:       if (ind[0] == max-1) {ind[0] = -1;}
173:       else                 {++ind[0];}
174:     }
175:   }
176:   return(0);
177: }

179: /*
180:   p in [0, npoints), i in [0, pdim), c in [0, Nc)

182:   B[p][i][c] = B[p][i_scalar][c][c]
183: */
184: PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
185: {
186:   PetscSpace_Poly *poly    = (PetscSpace_Poly *) sp->data;
187:   DM               dm      = sp->dm;
188:   PetscInt         Nc      = sp->Nc;
189:   PetscInt         ndegree = sp->degree+1;
190:   PetscInt        *degrees = poly->degrees;
191:   PetscInt         dim     = sp->Nv;
192:   PetscReal       *lpoints, *tmp, *LB, *LD, *LH;
193:   PetscInt        *ind, *tup;
194:   PetscInt         c, pdim, d, e, der, der2, i, p, deg, o;
195:   PetscErrorCode   ierr;

198:   PetscSpaceGetDimension(sp, &pdim);
199:   pdim /= Nc;
200:   DMGetWorkArray(dm, npoints, MPIU_REAL, &lpoints);
201:   DMGetWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);
202:   if (B || D || H) {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);}
203:   if (D || H)      {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);}
204:   if (H)           {DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);}
205:   for (d = 0; d < dim; ++d) {
206:     for (p = 0; p < npoints; ++p) {
207:       lpoints[p] = points[p*dim+d];
208:     }
209:     PetscDTLegendreEval(npoints, lpoints, ndegree, degrees, tmp, &tmp[1*npoints*ndegree], &tmp[2*npoints*ndegree]);
210:     /* LB, LD, LH (ndegree * dim x npoints) */
211:     for (deg = 0; deg < ndegree; ++deg) {
212:       for (p = 0; p < npoints; ++p) {
213:         if (B || D || H) LB[(deg*dim + d)*npoints + p] = tmp[(0*npoints + p)*ndegree+deg];
214:         if (D || H)      LD[(deg*dim + d)*npoints + p] = tmp[(1*npoints + p)*ndegree+deg];
215:         if (H)           LH[(deg*dim + d)*npoints + p] = tmp[(2*npoints + p)*ndegree+deg];
216:       }
217:     }
218:   }
219:   /* Multiply by A (pdim x ndegree * dim) */
220:   PetscMalloc2(dim,&ind,dim,&tup);
221:   if (B) {
222:     /* B (npoints x pdim x Nc) */
223:     PetscMemzero(B, npoints*pdim*Nc*Nc * sizeof(PetscReal));
224:     if (poly->tensor) {
225:       i = 0;
226:       PetscMemzero(ind, dim * sizeof(PetscInt));
227:       while (ind[0] >= 0) {
228:         TensorPoint_Internal(dim, sp->degree+1, ind, tup);
229:         for (p = 0; p < npoints; ++p) {
230:           B[(p*pdim + i)*Nc*Nc] = 1.0;
231:           for (d = 0; d < dim; ++d) {
232:             B[(p*pdim + i)*Nc*Nc] *= LB[(tup[d]*dim + d)*npoints + p];
233:           }
234:         }
235:         ++i;
236:       }
237:     } else {
238:       i = 0;
239:       for (o = 0; o <= sp->degree; ++o) {
240:         PetscMemzero(ind, dim * sizeof(PetscInt));
241:         while (ind[0] >= 0) {
242:           LatticePoint_Internal(dim, o, ind, tup);
243:           for (p = 0; p < npoints; ++p) {
244:             B[(p*pdim + i)*Nc*Nc] = 1.0;
245:             for (d = 0; d < dim; ++d) {
246:               B[(p*pdim + i)*Nc*Nc] *= LB[(tup[d]*dim + d)*npoints + p];
247:             }
248:           }
249:           ++i;
250:         }
251:       }
252:     }
253:     /* Make direct sum basis for multicomponent space */
254:     for (p = 0; p < npoints; ++p) {
255:       for (i = 0; i < pdim; ++i) {
256:         for (c = 1; c < Nc; ++c) {
257:           B[(p*pdim*Nc + i*Nc + c)*Nc + c] = B[(p*pdim + i)*Nc*Nc];
258:         }
259:       }
260:     }
261:   }
262:   if (D) {
263:     /* D (npoints x pdim x Nc x dim) */
264:     PetscMemzero(D, npoints*pdim*Nc*Nc*dim * sizeof(PetscReal));
265:     if (poly->tensor) {
266:       i = 0;
267:       PetscMemzero(ind, dim * sizeof(PetscInt));
268:       while (ind[0] >= 0) {
269:         TensorPoint_Internal(dim, sp->degree+1, ind, tup);
270:         for (p = 0; p < npoints; ++p) {
271:           for (der = 0; der < dim; ++der) {
272:             D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0;
273:             for (d = 0; d < dim; ++d) {
274:               if (d == der) {
275:                 D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
276:               } else {
277:                 D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
278:               }
279:             }
280:           }
281:         }
282:         ++i;
283:       }
284:     } else {
285:       i = 0;
286:       for (o = 0; o <= sp->degree; ++o) {
287:         PetscMemzero(ind, dim * sizeof(PetscInt));
288:         while (ind[0] >= 0) {
289:           LatticePoint_Internal(dim, o, ind, tup);
290:           for (p = 0; p < npoints; ++p) {
291:             for (der = 0; der < dim; ++der) {
292:               D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0;
293:               for (d = 0; d < dim; ++d) {
294:                 if (d == der) {
295:                   D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
296:                 } else {
297:                   D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
298:                 }
299:               }
300:             }
301:           }
302:           ++i;
303:         }
304:       }
305:     }
306:     /* Make direct sum basis for multicomponent space */
307:     for (p = 0; p < npoints; ++p) {
308:       for (i = 0; i < pdim; ++i) {
309:         for (c = 1; c < Nc; ++c) {
310:           for (d = 0; d < dim; ++d) {
311:             D[((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d] = D[(p*pdim + i)*Nc*Nc*dim + d];
312:           }
313:         }
314:       }
315:     }
316:   }
317:   if (H) {
318:     /* H (npoints x pdim x Nc x Nc x dim x dim) */
319:     PetscMemzero(H, npoints*pdim*Nc*Nc*dim*dim * sizeof(PetscReal));
320:     if (poly->tensor) {
321:       i = 0;
322:       PetscMemzero(ind, dim * sizeof(PetscInt));
323:       while (ind[0] >= 0) {
324:         TensorPoint_Internal(dim, sp->degree+1, ind, tup);
325:         for (p = 0; p < npoints; ++p) {
326:           for (der = 0; der < dim; ++der) {
327:             H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] = 1.0;
328:             for (d = 0; d < dim; ++d) {
329:               if (d == der) {
330:                 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LH[(tup[d]*dim + d)*npoints + p];
331:               } else {
332:                 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
333:               }
334:             }
335:             for (der2 = der + 1; der2 < dim; ++der2) {
336:               H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0;
337:               for (d = 0; d < dim; ++d) {
338:                 if (d == der || d == der2) {
339:                   H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p];
340:                 } else {
341:                   H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p];
342:                 }
343:               }
344:               H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2];
345:             }
346:           }
347:         }
348:         ++i;
349:       }
350:     } else {
351:       i = 0;
352:       for (o = 0; o <= sp->degree; ++o) {
353:         PetscMemzero(ind, dim * sizeof(PetscInt));
354:         while (ind[0] >= 0) {
355:           LatticePoint_Internal(dim, o, ind, tup);
356:           for (p = 0; p < npoints; ++p) {
357:             for (der = 0; der < dim; ++der) {
358:               H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] = 1.0;
359:               for (d = 0; d < dim; ++d) {
360:                 if (d == der) {
361:                   H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LH[(tup[d]*dim + d)*npoints + p];
362:                 } else {
363:                   H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
364:                 }
365:               }
366:               for (der2 = der + 1; der2 < dim; ++der2) {
367:                 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0;
368:                 for (d = 0; d < dim; ++d) {
369:                   if (d == der || d == der2) {
370:                     H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p];
371:                   } else {
372:                     H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p];
373:                   }
374:                 }
375:                 H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2];
376:               }
377:             }
378:           }
379:           ++i;
380:         }
381:       }
382:     }
383:     /* Make direct sum basis for multicomponent space */
384:     for (p = 0; p < npoints; ++p) {
385:       for (i = 0; i < pdim; ++i) {
386:         for (c = 1; c < Nc; ++c) {
387:           for (d = 0; d < dim; ++d) {
388:             for (e = 0; e < dim; ++e) {
389:               H[(((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d)*dim + e] = H[((p*pdim + i)*Nc*Nc*dim + d)*dim + e];
390:             }
391:           }
392:         }
393:       }
394:     }
395:   }
396:   PetscFree2(ind,tup);
397:   if (H)           {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);}
398:   if (D || H)      {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);}
399:   if (B || D || H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);}
400:   DMRestoreWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);
401:   DMRestoreWorkArray(dm, npoints, MPIU_REAL, &lpoints);
402:   return(0);
403: }

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

410:   Input Parameters:
411: + sp     - the function space object
412: - tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space

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

417:   Level: beginner

419: .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
420: @*/
421: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
422: {

427:   PetscTryMethod(sp,"PetscSpacePolynomialSetTensor_C",(PetscSpace,PetscBool),(sp,tensor));
428:   return(0);
429: }

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

436:   Input Parameters:
437: . sp     - the function space object

439:   Output Parameters:
440: . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space

442:   Level: beginner

444: .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
445: @*/
446: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
447: {

453:   PetscTryMethod(sp,"PetscSpacePolynomialGetTensor_C",(PetscSpace,PetscBool*),(sp,tensor));
454:   return(0);
455: }

457: static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
458: {
459:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

462:   poly->tensor = tensor;
463:   return(0);
464: }

466: static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor)
467: {
468:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

473:   *tensor = poly->tensor;
474:   return(0);
475: }

477: static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp)
478: {
479:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
480:   PetscInt         Nc, dim, order;
481:   PetscBool        tensor;
482:   PetscErrorCode   ierr;

485:   PetscSpaceGetNumComponents(sp, &Nc);
486:   PetscSpaceGetNumVariables(sp, &dim);
487:   PetscSpaceGetDegree(sp, &order, NULL);
488:   PetscSpacePolynomialGetTensor(sp, &tensor);
489:   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);}
490:   if (!poly->subspaces) {PetscCalloc1(dim, &poly->subspaces);}
491:   if (height <= dim) {
492:     if (!poly->subspaces[height-1]) {
493:       PetscSpace  sub;
494:       const char *name;

496:       PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
497:       PetscObjectGetName((PetscObject) sp,  &name);
498:       PetscObjectSetName((PetscObject) sub,  name);
499:       PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL);
500:       PetscSpaceSetNumComponents(sub, Nc);
501:       PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
502:       PetscSpaceSetNumVariables(sub, dim-height);
503:       PetscSpacePolynomialSetTensor(sub, tensor);
504:       PetscSpaceSetUp(sub);
505:       poly->subspaces[height-1] = sub;
506:     }
507:     *subsp = poly->subspaces[height-1];
508:   } else {
509:     *subsp = NULL;
510:   }
511:   return(0);
512: }

514: PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
515: {

519:   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_Polynomial;
520:   sp->ops->setup             = PetscSpaceSetUp_Polynomial;
521:   sp->ops->view              = PetscSpaceView_Polynomial;
522:   sp->ops->destroy           = PetscSpaceDestroy_Polynomial;
523:   sp->ops->getdimension      = PetscSpaceGetDimension_Polynomial;
524:   sp->ops->evaluate          = PetscSpaceEvaluate_Polynomial;
525:   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial;
526:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial);
527:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial);
528:   return(0);
529: }

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

535:   Level: intermediate

537: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
538: M*/

540: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
541: {
542:   PetscSpace_Poly *poly;
543:   PetscErrorCode   ierr;

547:   PetscNewLog(sp,&poly);
548:   sp->data = poly;

550:   poly->symmetric    = PETSC_FALSE;
551:   poly->tensor       = PETSC_FALSE;
552:   poly->degrees      = NULL;
553:   poly->subspaces    = NULL;

555:   PetscSpaceInitialize_Polynomial(sp);
556:   return(0);
557: }

559: PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym)
560: {
561:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

565:   poly->symmetric = sym;
566:   return(0);
567: }

569: PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym)
570: {
571:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

576:   *sym = poly->symmetric;
577:   return(0);
578: }