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, pStart, pEnd;
  9:   PetscSection           section;

 11:   DMGetDimension(dm, &dim);
 12:   DMPlexGetChart(dm, &pStart, &pEnd);
 13:   PetscSectionCreate(PETSC_COMM_SELF, &section);
 14:   PetscSectionSetChart(section, pStart, pEnd);
 15:   PetscSectionSetDof(section, pStart, s->dim);
 16:   PetscSectionSetUp(section);
 17:   sp->pointSection = section;
 18:   return 0;
 19: }

 21: static PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
 22: {
 23:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;

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

 32: static PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace spNew)
 33: {
 34:   PetscInt       dim, d;

 36:   PetscDualSpaceGetDimension(sp, &dim);
 37:   PetscDualSpaceSimpleSetDimension(spNew, dim);
 38:   for (d = 0; d < dim; ++d) {
 39:     PetscQuadrature q;

 41:     PetscDualSpaceGetFunctional(sp, d, &q);
 42:     PetscDualSpaceSimpleSetFunctional(spNew, d, q);
 43:   }
 44:   return 0;
 45: }

 47: static PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
 48: {
 49:   return 0;
 50: }

 52: static PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
 53: {
 54:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
 55:   DM                     dm;
 56:   PetscInt               spatialDim, f;

 58:   for (f = 0; f < s->dim; ++f) PetscQuadratureDestroy(&sp->functional[f]);
 59:   PetscFree(sp->functional);
 60:   s->dim = dim;
 61:   PetscCalloc1(s->dim, &sp->functional);
 62:   PetscFree(s->numDof);
 63:   PetscDualSpaceGetDM(sp, &dm);
 64:   DMGetCoordinateDim(dm, &spatialDim);
 65:   PetscCalloc1(spatialDim+1, &s->numDof);
 66:   s->numDof[spatialDim] = dim;
 67:   return 0;
 68: }

 70: static PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
 71: {
 72:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
 73:   PetscReal             *weights;
 74:   PetscInt               Nc, c, Nq, p;

 77:   PetscQuadratureDuplicate(q, &sp->functional[f]);
 78:   /* Reweight so that it has unit volume: Do we want to do this for Nc > 1? */
 79:   PetscQuadratureGetData(sp->functional[f], NULL, &Nc, &Nq, NULL, (const PetscReal **) &weights);
 80:   for (c = 0; c < Nc; ++c) {
 81:     PetscReal vol = 0.0;

 83:     for (p = 0; p < Nq; ++p) vol += weights[p*Nc+c];
 84:     for (p = 0; p < Nq; ++p) weights[p*Nc+c] /= (vol == 0.0 ? 1.0 : vol);
 85:   }
 86:   return 0;
 87: }

 89: /*@
 90:   PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis

 92:   Logically Collective on sp

 94:   Input Parameters:
 95: + sp  - the PetscDualSpace
 96: - dim - the basis dimension

 98:   Level: intermediate

100: .seealso: PetscDualSpaceSimpleSetFunctional()
101: @*/
102: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
103: {
107:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));
108:   return 0;
109: }

111: /*@
112:   PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space

114:   Not Collective

116:   Input Parameters:
117: + sp  - the PetscDualSpace
118: . f - the basis index
119: - q - the basis functional

121:   Level: intermediate

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

125: .seealso: PetscDualSpaceSimpleSetDimension()
126: @*/
127: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
128: {
130:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));
131:   return 0;
132: }

134: static PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
135: {
136:   sp->ops->setfromoptions       = PetscDualSpaceSetFromOptions_Simple;
137:   sp->ops->setup                = PetscDualSpaceSetUp_Simple;
138:   sp->ops->view                 = NULL;
139:   sp->ops->destroy              = PetscDualSpaceDestroy_Simple;
140:   sp->ops->duplicate            = PetscDualSpaceDuplicate_Simple;
141:   sp->ops->createheightsubspace = NULL;
142:   sp->ops->createpointsubspace  = NULL;
143:   sp->ops->getsymmetries        = NULL;
144:   sp->ops->apply                = PetscDualSpaceApplyDefault;
145:   sp->ops->applyall             = PetscDualSpaceApplyAllDefault;
146:   sp->ops->applyint             = PetscDualSpaceApplyInteriorDefault;
147:   sp->ops->createalldata        = PetscDualSpaceCreateAllDataDefault;
148:   sp->ops->createintdata        = PetscDualSpaceCreateInteriorDataDefault;
149:   return 0;
150: }

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

155:   Level: intermediate

157: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
158: M*/

160: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
161: {
162:   PetscDualSpace_Simple *s;

165:   PetscNewLog(sp,&s);
166:   sp->data = s;

168:   s->dim    = 0;
169:   s->numDof = NULL;

171:   PetscDualSpaceInitialize_Simple(sp);
172:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);
173:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);
174:   return 0;
175: }