1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petscdmplex.h> 3 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 11 PetscFunctionBegin; 12 PetscCall(DMGetDimension(dm, &dim)); 13 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 14 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, §ion)); 15 PetscCall(PetscSectionSetChart(section, pStart, pEnd)); 16 PetscCall(PetscSectionSetDof(section, pStart, s->dim)); 17 PetscCall(PetscSectionSetUp(section)); 18 sp->pointSection = section; 19 PetscFunctionReturn(0); 20 } 21 22 static PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp) 23 { 24 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data; 25 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(0); 32 } 33 34 static PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace spNew) 35 { 36 PetscInt dim, d; 37 38 PetscFunctionBegin; 39 PetscCall(PetscDualSpaceGetDimension(sp, &dim)); 40 PetscCall(PetscDualSpaceSimpleSetDimension(spNew, dim)); 41 for (d = 0; d < dim; ++d) { 42 PetscQuadrature q; 43 44 PetscCall(PetscDualSpaceGetFunctional(sp, d, &q)); 45 PetscCall(PetscDualSpaceSimpleSetFunctional(spNew, d, q)); 46 } 47 PetscFunctionReturn(0); 48 } 49 50 static PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp) 51 { 52 PetscFunctionBegin; 53 PetscFunctionReturn(0); 54 } 55 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; 61 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(0); 73 } 74 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; 80 81 PetscFunctionBegin; 82 PetscCheckFalse((f < 0) || (f >= s->dim),PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %d not in [0, %d)", 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; 88 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(0); 93 } 94 95 /*@ 96 PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis 97 98 Logically Collective on sp 99 100 Input Parameters: 101 + sp - the PetscDualSpace 102 - dim - the basis dimension 103 104 Level: intermediate 105 106 .seealso: PetscDualSpaceSimpleSetFunctional() 107 @*/ 108 PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim) 109 { 110 PetscFunctionBegin; 111 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 112 PetscValidLogicalCollectiveInt(sp, dim, 2); 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(0); 116 } 117 118 /*@ 119 PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space 120 121 Not Collective 122 123 Input Parameters: 124 + sp - the PetscDualSpace 125 . f - the basis index 126 - q - the basis functional 127 128 Level: intermediate 129 130 Note: The quadrature will be reweighted so that it has unit volume. 131 132 .seealso: PetscDualSpaceSimpleSetDimension() 133 @*/ 134 PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q) 135 { 136 PetscFunctionBegin; 137 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 138 PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q)); 139 PetscFunctionReturn(0); 140 } 141 142 static PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp) 143 { 144 PetscFunctionBegin; 145 sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Simple; 146 sp->ops->setup = PetscDualSpaceSetUp_Simple; 147 sp->ops->view = NULL; 148 sp->ops->destroy = PetscDualSpaceDestroy_Simple; 149 sp->ops->duplicate = PetscDualSpaceDuplicate_Simple; 150 sp->ops->createheightsubspace = NULL; 151 sp->ops->createpointsubspace = NULL; 152 sp->ops->getsymmetries = NULL; 153 sp->ops->apply = PetscDualSpaceApplyDefault; 154 sp->ops->applyall = PetscDualSpaceApplyAllDefault; 155 sp->ops->applyint = PetscDualSpaceApplyInteriorDefault; 156 sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault; 157 sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault; 158 PetscFunctionReturn(0); 159 } 160 161 /*MC 162 PETSCDUALSPACESIMPLE = "simple" - A PetscDualSpace object that encapsulates a dual space of arbitrary functionals 163 164 Level: intermediate 165 166 .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType() 167 M*/ 168 169 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp) 170 { 171 PetscDualSpace_Simple *s; 172 173 PetscFunctionBegin; 174 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 175 PetscCall(PetscNewLog(sp,&s)); 176 sp->data = s; 177 178 s->dim = 0; 179 s->numDof = NULL; 180 181 PetscCall(PetscDualSpaceInitialize_Simple(sp)); 182 PetscCall(PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple)); 183 PetscCall(PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple)); 184 PetscFunctionReturn(0); 185 } 186