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

  4: static 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: static 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: static 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: static PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
 52: {
 54:   return(0);
 55: }

 57: static 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: static 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: static 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: static 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 sp

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

125:   Level: intermediate

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

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

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

143:   Not Collective

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

150:   Level: intermediate

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

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

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

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

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

187:   Level: intermediate

189: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
190: M*/

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

199:   PetscNewLog(sp,&s);
200:   sp->data = s;

202:   s->dim    = 0;
203:   s->numDof = NULL;

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