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;
 10:   PetscErrorCode         ierr;

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

 23: static PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
 24: {
 25:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
 26:   PetscErrorCode         ierr;

 29:   PetscFree(s->numDof);
 30:   PetscFree(s);
 31:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", NULL);
 32:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", NULL);
 33:   return(0);
 34: }

 36: static PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace spNew)
 37: {
 38:   PetscInt       dim, d;

 42:   PetscDualSpaceGetDimension(sp, &dim);
 43:   PetscDualSpaceSimpleSetDimension(spNew, dim);
 44:   for (d = 0; d < dim; ++d) {
 45:     PetscQuadrature q;

 47:     PetscDualSpaceGetFunctional(sp, d, &q);
 48:     PetscDualSpaceSimpleSetFunctional(spNew, d, q);
 49:   }
 50:   return(0);
 51: }

 53: static PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
 54: {
 56:   return(0);
 57: }

 59: static PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
 60: {
 61:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
 62:   DM                     dm;
 63:   PetscInt               spatialDim, f;
 64:   PetscErrorCode         ierr;

 67:   for (f = 0; f < s->dim; ++f) {PetscQuadratureDestroy(&sp->functional[f]);}
 68:   PetscFree(sp->functional);
 69:   s->dim = dim;
 70:   PetscCalloc1(s->dim, &sp->functional);
 71:   PetscFree(s->numDof);
 72:   PetscDualSpaceGetDM(sp, &dm);
 73:   DMGetCoordinateDim(dm, &spatialDim);
 74:   PetscCalloc1(spatialDim+1, &s->numDof);
 75:   s->numDof[spatialDim] = dim;
 76:   return(0);
 77: }

 79: static PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
 80: {
 81:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
 82:   PetscReal             *weights;
 83:   PetscInt               Nc, c, Nq, p;
 84:   PetscErrorCode         ierr;

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

 94:     for (p = 0; p < Nq; ++p) vol += weights[p*Nc+c];
 95:     for (p = 0; p < Nq; ++p) weights[p*Nc+c] /= (vol == 0.0 ? 1.0 : vol);
 96:   }
 97:   return(0);
 98: }

100: /*@
101:   PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis

103:   Logically Collective on sp

105:   Input Parameters:
106: + sp  - the PetscDualSpace
107: - dim - the basis dimension

109:   Level: intermediate

111: .seealso: PetscDualSpaceSimpleSetFunctional()
112: @*/
113: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
114: {

120:   if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change dimension after dual space is set up");
121:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));
122:   return(0);
123: }

125: /*@
126:   PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space

128:   Not Collective

130:   Input Parameters:
131: + sp  - the PetscDualSpace
132: . f - the basis index
133: - q - the basis functional

135:   Level: intermediate

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

139: .seealso: PetscDualSpaceSimpleSetDimension()
140: @*/
141: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
142: {

147:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));
148:   return(0);
149: }

151: static PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
152: {
154:   sp->ops->setfromoptions       = PetscDualSpaceSetFromOptions_Simple;
155:   sp->ops->setup                = PetscDualSpaceSetUp_Simple;
156:   sp->ops->view                 = NULL;
157:   sp->ops->destroy              = PetscDualSpaceDestroy_Simple;
158:   sp->ops->duplicate            = PetscDualSpaceDuplicate_Simple;
159:   sp->ops->createheightsubspace = NULL;
160:   sp->ops->createpointsubspace  = NULL;
161:   sp->ops->getsymmetries        = NULL;
162:   sp->ops->apply                = PetscDualSpaceApplyDefault;
163:   sp->ops->applyall             = PetscDualSpaceApplyAllDefault;
164:   sp->ops->applyint             = PetscDualSpaceApplyInteriorDefault;
165:   sp->ops->createalldata        = PetscDualSpaceCreateAllDataDefault;
166:   sp->ops->createintdata        = PetscDualSpaceCreateInteriorDataDefault;
167:   return(0);
168: }

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

173:   Level: intermediate

175: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
176: M*/

178: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
179: {
180:   PetscDualSpace_Simple *s;
181:   PetscErrorCode         ierr;

185:   PetscNewLog(sp,&s);
186:   sp->data = s;

188:   s->dim    = 0;
189:   s->numDof = NULL;

191:   PetscDualSpaceInitialize_Simple(sp);
192:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);
193:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);
194:   return(0);
195: }