Actual source code: dspacesimple.c

  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(PetscDualSpace sp, PetscOptionItems *PetscOptionsObject)
 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: `PetscDualSpace`, `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:
124:   The quadrature will be reweighted so that it has unit volume.

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

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

153: /*MC
154:   PETSCDUALSPACESIMPLE = "simple" - A `PetscDualSpaceType` that encapsulates a dual space of arbitrary functionals

156:   Level: intermediate

158: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
159: M*/

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

166:   PetscNew(&s);
167:   sp->data = s;

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

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