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 PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
 51: {
 52:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data;
 53:   DM                     dm;
 54:   PetscInt               spatialDim, f;

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

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

 75:   PetscFunctionBegin;
 76:   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);
 77:   PetscCall(PetscQuadratureDuplicate(q, &sp->functional[f]));
 78:   /* Reweight so that it has unit volume: Do we want to do this for Nc > 1? */
 79:   PetscCall(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:   PetscFunctionReturn(PETSC_SUCCESS);
 87: }

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

 92:   Logically Collective

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

 98:   Level: intermediate

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

112: /*@
113:   PetscDualSpaceSimpleSetFunctional - Set the given basis functional for this dual space

115:   Not Collective

117:   Input Parameters:
118: + sp   - the `PetscDualSpace`
119: . func - the basis index
120: - q    - the basis functional

122:   Level: intermediate

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

127: .seealso: `PETSCDUALSPACESIMPLE`, `PetscDualSpace`, `PetscDualSpaceSimpleSetDimension()`
128: @*/
129: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
130: {
131:   PetscFunctionBegin;
133:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace, PetscInt, PetscQuadrature), (sp, func, q));
134:   PetscFunctionReturn(PETSC_SUCCESS);
135: }

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

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

158:   Level: intermediate

160:   Developer Note:
161:   It is not clear this has a good name

163: .seealso: `PetscDualSpace`, `PetscDualSpaceSimpleSetFunctional()`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
164: M*/

166: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
167: {
168:   PetscDualSpace_Simple *s;

170:   PetscFunctionBegin;
172:   PetscCall(PetscNew(&s));
173:   sp->data = s;

175:   s->dim    = 0;
176:   s->numDof = NULL;

178:   PetscCall(PetscDualSpaceInitialize_Simple(sp));
179:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple));
180:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple));
181:   PetscFunctionReturn(PETSC_SUCCESS);
182: }