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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 48 } 49 50 static PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscDualSpace sp, PetscOptionItems *PetscOptionsObject) 51 { 52 PetscFunctionBegin; 53 PetscFunctionReturn(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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 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; 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(PETSC_SUCCESS); 93 } 94 95 /*@ 96 PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis 97 98 Logically Collective 99 100 Input Parameters: 101 + sp - the `PetscDualSpace` 102 - dim - the basis dimension 103 104 Level: intermediate 105 106 .seealso: `PETSCDUALSPACESIMPLE`, `PetscDualSpace`, `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(PETSC_SUCCESS); 116 } 117 118 /*@ 119 PetscDualSpaceSimpleSetFunctional - Set the given basis functional for this dual space 120 121 Not Collective 122 123 Input Parameters: 124 + sp - the `PetscDualSpace` 125 . func - the basis index 126 - q - the basis functional 127 128 Level: intermediate 129 130 Note: 131 The quadrature will be reweighted so that it has unit volume. 132 133 .seealso: `PETSCDUALSPACESIMPLE`, `PetscDualSpace`, `PetscDualSpaceSimpleSetDimension()` 134 @*/ 135 PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q) 136 { 137 PetscFunctionBegin; 138 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 139 PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace, PetscInt, PetscQuadrature), (sp, func, q)); 140 PetscFunctionReturn(PETSC_SUCCESS); 141 } 142 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 } 161 162 /*MC 163 PETSCDUALSPACESIMPLE = "simple" - A `PetscDualSpaceType` that encapsulates a dual space of functionals provided with `PetscDualSpaceSimpleSetFunctional()` 164 165 Level: intermediate 166 167 Developer Note: 168 It is not clear this has a good name 169 170 .seealso: `PetscDualSpace`, `PetscDualSpaceSimpleSetFunctional()`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()` 171 M*/ 172 173 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp) 174 { 175 PetscDualSpace_Simple *s; 176 177 PetscFunctionBegin; 178 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 179 PetscCall(PetscNew(&s)); 180 sp->data = s; 181 182 s->dim = 0; 183 s->numDof = NULL; 184 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 } 190