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:   PetscFunctionBegin;
 12:   PetscCall(DMGetDimension(dm, &dim));
 13:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
 14:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &section));
 15:   PetscCall(PetscSectionSetChart(section, pStart, pEnd));
 16:   PetscCall(PetscSectionSetDof(section, pStart, s->dim));
 17:   PetscCall(PetscSectionSetUp(section));
 18:   sp->pointSection = section;
 19:   PetscFunctionReturn(PETSC_SUCCESS);
 20: }

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

 26:   PetscFunctionBegin;
 27:   PetscCall(PetscFree(s->numDof));
 28:   PetscCall(PetscFree(s));
 29:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetDimension_C", NULL));
 30:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetFunctional_C", NULL));
 31:   PetscFunctionReturn(PETSC_SUCCESS);
 32: }

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

 38:   PetscFunctionBegin;
 39:   PetscCall(PetscDualSpaceGetDimension(sp, &dim));
 40:   PetscCall(PetscDualSpaceSimpleSetDimension(spNew, dim));
 41:   for (d = 0; d < dim; ++d) {
 42:     PetscQuadrature q;

 44:     PetscCall(PetscDualSpaceGetFunctional(sp, d, &q));
 45:     PetscCall(PetscDualSpaceSimpleSetFunctional(spNew, d, q));
 46:   }
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

 50: static PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscDualSpace sp, PetscOptionItems *PetscOptionsObject)
 51: {
 52:   PetscFunctionBegin;
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

 56: static PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
 57: {
 58:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
 59:   DM                     dm;
 60:   PetscInt               spatialDim, f;

 62:   PetscFunctionBegin;
 63:   for (f = 0; f < s->dim; ++f) PetscCall(PetscQuadratureDestroy(&sp->functional[f]));
 64:   PetscCall(PetscFree(sp->functional));
 65:   s->dim = dim;
 66:   PetscCall(PetscCalloc1(s->dim, &sp->functional));
 67:   PetscCall(PetscFree(s->numDof));
 68:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
 69:   PetscCall(DMGetCoordinateDim(dm, &spatialDim));
 70:   PetscCall(PetscCalloc1(spatialDim + 1, &s->numDof));
 71:   s->numDof[spatialDim] = dim;
 72:   PetscFunctionReturn(PETSC_SUCCESS);
 73: }

 75: static PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
 76: {
 77:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
 78:   PetscReal             *weights;
 79:   PetscInt               Nc, c, Nq, p;

 81:   PetscFunctionBegin;
 82:   PetscCheck(!(f < 0) && !(f >= s->dim), PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", f, s->dim);
 83:   PetscCall(PetscQuadratureDuplicate(q, &sp->functional[f]));
 84:   /* Reweight so that it has unit volume: Do we want to do this for Nc > 1? */
 85:   PetscCall(PetscQuadratureGetData(sp->functional[f], NULL, &Nc, &Nq, NULL, (const PetscReal **)&weights));
 86:   for (c = 0; c < Nc; ++c) {
 87:     PetscReal vol = 0.0;

 89:     for (p = 0; p < Nq; ++p) vol += weights[p * Nc + c];
 90:     for (p = 0; p < Nq; ++p) weights[p * Nc + c] /= (vol == 0.0 ? 1.0 : vol);
 91:   }
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 93: }

 95: /*@
 96:   PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis

 98:   Logically Collective

100:   Input Parameters:
101: + sp  - the `PetscDualSpace`
102: - dim - the basis dimension

104:   Level: intermediate

106: .seealso: `PETSCDUALSPACESIMPLE`, `PetscDualSpace`, `PetscDualSpaceSimpleSetFunctional()`
107: @*/
108: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
109: {
110:   PetscFunctionBegin;
113:   PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change dimension after dual space is set up");
114:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace, PetscInt), (sp, dim));
115:   PetscFunctionReturn(PETSC_SUCCESS);
116: }

118: /*@
119:   PetscDualSpaceSimpleSetFunctional - Set the given basis functional for this dual space

121:   Not Collective

123:   Input Parameters:
124: + sp   - the `PetscDualSpace`
125: . func - the basis index
126: - q    - the basis functional

128:   Level: intermediate

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

133: .seealso: `PETSCDUALSPACESIMPLE`, `PetscDualSpace`, `PetscDualSpaceSimpleSetDimension()`
134: @*/
135: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
136: {
137:   PetscFunctionBegin;
139:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace, PetscInt, PetscQuadrature), (sp, func, q));
140:   PetscFunctionReturn(PETSC_SUCCESS);
141: }

143: static PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
144: {
145:   PetscFunctionBegin;
146:   sp->ops->setfromoptions       = PetscDualSpaceSetFromOptions_Simple;
147:   sp->ops->setup                = PetscDualSpaceSetUp_Simple;
148:   sp->ops->view                 = NULL;
149:   sp->ops->destroy              = PetscDualSpaceDestroy_Simple;
150:   sp->ops->duplicate            = PetscDualSpaceDuplicate_Simple;
151:   sp->ops->createheightsubspace = NULL;
152:   sp->ops->createpointsubspace  = NULL;
153:   sp->ops->getsymmetries        = NULL;
154:   sp->ops->apply                = PetscDualSpaceApplyDefault;
155:   sp->ops->applyall             = PetscDualSpaceApplyAllDefault;
156:   sp->ops->applyint             = PetscDualSpaceApplyInteriorDefault;
157:   sp->ops->createalldata        = PetscDualSpaceCreateAllDataDefault;
158:   sp->ops->createintdata        = PetscDualSpaceCreateInteriorDataDefault;
159:   PetscFunctionReturn(PETSC_SUCCESS);
160: }

162: /*MC
163:   PETSCDUALSPACESIMPLE = "simple" - A `PetscDualSpaceType` that encapsulates a dual space of functionals provided with `PetscDualSpaceSimpleSetFunctional()`

165:   Level: intermediate

167:   Developer Note:
168:   It is not clear this has a good name

170: .seealso: `PetscDualSpace`, `PetscDualSpaceSimpleSetFunctional()`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
171: M*/

173: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
174: {
175:   PetscDualSpace_Simple *s;

177:   PetscFunctionBegin;
179:   PetscCall(PetscNew(&s));
180:   sp->data = s;

182:   s->dim    = 0;
183:   s->numDof = NULL;

185:   PetscCall(PetscDualSpaceInitialize_Simple(sp));
186:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple));
187:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple));
188:   PetscFunctionReturn(PETSC_SUCCESS);
189: }