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

  4: PetscErrorCode PetscDualSpaceSetUp_Simple(PetscDualSpace sp)
  5: {
  6:   PetscDualSpace_Simple *s  = (PetscDualSpace_Simple *) sp->data;
  7:   DM                     dm = sp->dm;
  8:   PetscInt               dim;
  9:   PetscErrorCode         ierr;

 12:   DMGetDimension(dm, &dim);
 13:   PetscCalloc1(dim+1, &s->numDof);
 14:   return(0);
 15: }

 17: PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
 18: {
 19:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
 20:   PetscErrorCode         ierr;

 23:   PetscFree(s->numDof);
 24:   PetscFree(s);
 25:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", NULL);
 26:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", NULL);
 27:   return(0);
 28: }

 30: PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace *spNew)
 31: {
 32:   PetscInt       dim, d, Nc;

 36:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
 37:   PetscDualSpaceSetType(*spNew, PETSCDUALSPACESIMPLE);
 38:   PetscDualSpaceGetNumComponents(sp, &Nc);
 39:   PetscDualSpaceSetNumComponents(sp, Nc);
 40:   PetscDualSpaceGetDimension(sp, &dim);
 41:   PetscDualSpaceSimpleSetDimension(*spNew, dim);
 42:   for (d = 0; d < dim; ++d) {
 43:     PetscQuadrature q;

 45:     PetscDualSpaceGetFunctional(sp, d, &q);
 46:     PetscDualSpaceSimpleSetFunctional(*spNew, d, q);
 47:   }
 48:   return(0);
 49: }

 51: PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
 52: {
 54:   return(0);
 55: }

 57: PetscErrorCode PetscDualSpaceGetDimension_Simple(PetscDualSpace sp, PetscInt *dim)
 58: {
 59:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;

 62:   *dim = s->dim;
 63:   return(0);
 64: }

 66: PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
 67: {
 68:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
 69:   DM                     dm;
 70:   PetscInt               spatialDim, f;
 71:   PetscErrorCode         ierr;

 74:   for (f = 0; f < s->dim; ++f) {PetscQuadratureDestroy(&sp->functional[f]);}
 75:   PetscFree(sp->functional);
 76:   s->dim = dim;
 77:   PetscCalloc1(s->dim, &sp->functional);
 78:   PetscFree(s->numDof);
 79:   PetscDualSpaceGetDM(sp, &dm);
 80:   DMGetCoordinateDim(dm, &spatialDim);
 81:   PetscCalloc1(spatialDim+1, &s->numDof);
 82:   s->numDof[spatialDim] = dim;
 83:   return(0);
 84: }

 86: PetscErrorCode PetscDualSpaceGetNumDof_Simple(PetscDualSpace sp, const PetscInt **numDof)
 87: {
 88:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;

 91:   *numDof = s->numDof;
 92:   return(0);
 93: }

 95: PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
 96: {
 97:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
 98:   PetscReal             *weights;
 99:   PetscInt               Nc, c, Nq, p;
100:   PetscErrorCode         ierr;

103:   if ((f < 0) || (f >= s->dim)) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %d not in [0, %d)", f, s->dim);
104:   PetscQuadratureDuplicate(q, &sp->functional[f]);
105:   /* Reweight so that it has unit volume: Do we want to do this for Nc > 1? */
106:   PetscQuadratureGetData(sp->functional[f], NULL, &Nc, &Nq, NULL, (const PetscReal **) &weights);
107:   for (c = 0; c < Nc; ++c) {
108:     PetscReal vol = 0.0;

110:     for (p = 0; p < Nq; ++p) vol += weights[p*Nc+c];
111:     for (p = 0; p < Nq; ++p) weights[p*Nc+c] /= (vol == 0.0 ? 1.0 : vol);
112:   }
113:   return(0);
114: }

116: /*@
117:   PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis

119:   Logically Collective on PetscDualSpace

121:   Input Parameters:
122: + sp  - the PetscDualSpace
123: - dim - the basis dimension

125:   Level: intermediate

127: .keywords: PetscDualSpace, dimension
128: .seealso: PetscDualSpaceSimpleSetFunctional()
129: @*/
130: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
131: {

137:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));
138:   return(0);
139: }

141: /*@
142:   PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space

144:   Not Collective

146:   Input Parameters:
147: + sp  - the PetscDualSpace
148: . f - the basis index
149: - q - the basis functional

151:   Level: intermediate

153:   Note: The quadrature will be reweighted so that it has unit volume.

155: .keywords: PetscDualSpace, functional
156: .seealso: PetscDualSpaceSimpleSetDimension()
157: @*/
158: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
159: {

164:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));
165:   return(0);
166: }

168: PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
169: {
171:   sp->ops->setfromoptions    = PetscDualSpaceSetFromOptions_Simple;
172:   sp->ops->setup             = PetscDualSpaceSetUp_Simple;
173:   sp->ops->view              = NULL;
174:   sp->ops->destroy           = PetscDualSpaceDestroy_Simple;
175:   sp->ops->duplicate         = PetscDualSpaceDuplicate_Simple;
176:   sp->ops->getdimension      = PetscDualSpaceGetDimension_Simple;
177:   sp->ops->getnumdof         = PetscDualSpaceGetNumDof_Simple;
178:   sp->ops->getheightsubspace = NULL;
179:   sp->ops->getsymmetries     = NULL;
180:   sp->ops->apply             = PetscDualSpaceApplyDefault;
181:   sp->ops->applyall          = PetscDualSpaceApplyAllDefault;
182:   sp->ops->createallpoints   = PetscDualSpaceCreateAllPointsDefault;
183:   return(0);
184: }

186: /*MC
187:   PETSCDUALSPACESIMPLE = "simple" - A PetscDualSpace object that encapsulates a dual space of arbitrary functionals

189:   Level: intermediate

191: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
192: M*/

194: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
195: {
196:   PetscDualSpace_Simple *s;
197:   PetscErrorCode         ierr;

201:   PetscNewLog(sp,&s);
202:   sp->data = s;

204:   s->dim    = 0;
205:   s->numDof = NULL;

207:   PetscDualSpaceInitialize_Simple(sp);
208:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);
209:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);
210:   return(0);
211: }