Actual source code: spacepoly.c

petsc-3.10.3 2018-12-18
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 viewer)
 17: {
 18:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
 19:   PetscErrorCode   ierr;

 22:   if (poly->tensor) {PetscViewerASCIIPrintf(viewer, "Tensor polynomial space of degree %D\n", sp->degree);}
 23:   else              {PetscViewerASCIIPrintf(viewer, "Polynomial space of degree %D\n", sp->degree);}
 24:   return(0);
 25: }

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

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

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

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

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

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

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

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

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

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

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

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

115:   Level: developer

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

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

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

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

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

153:   Level: developer

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

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

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

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

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

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

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

415:   Level: beginner

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

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

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

434:   Input Parameters:
435: . sp     - the function space object

437:   Output Parameters:
438: . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space

440:   Level: beginner

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

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

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

460:   poly->tensor = tensor;
461:   return(0);
462: }

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

471:   *tensor = poly->tensor;
472:   return(0);
473: }

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

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

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

509: PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
510: {

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

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

530:   Level: intermediate

532: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
533: M*/

535: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
536: {
537:   PetscSpace_Poly *poly;
538:   PetscErrorCode   ierr;

542:   PetscNewLog(sp,&poly);
543:   sp->data = poly;

545:   poly->symmetric    = PETSC_FALSE;
546:   poly->tensor       = PETSC_FALSE;
547:   poly->degrees      = NULL;
548:   poly->subspaces    = NULL;

550:   PetscSpaceInitialize_Polynomial(sp);
551:   return(0);
552: }

554: PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym)
555: {
556:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

560:   poly->symmetric = sym;
561:   return(0);
562: }

564: PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym)
565: {
566:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

571:   *sym = poly->symmetric;
572:   return(0);
573: }