Actual source code: dspacesimple.c
petsc-3.14.6 2021-03-30
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, §ion);
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: }