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 PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim) 51 { 52 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data; 53 DM dm; 54 PetscInt spatialDim, f; 55 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 } 68 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; 74 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; 82 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 } 88 89 /*@ 90 PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis 91 92 Logically Collective 93 94 Input Parameters: 95 + sp - the `PetscDualSpace` 96 - dim - the basis dimension 97 98 Level: intermediate 99 100 .seealso: `PETSCDUALSPACESIMPLE`, `PetscDualSpace`, `PetscDualSpaceSimpleSetFunctional()` 101 @*/ 102 PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim) 103 { 104 PetscFunctionBegin; 105 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 106 PetscValidLogicalCollectiveInt(sp, dim, 2); 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 } 111 112 /*@ 113 PetscDualSpaceSimpleSetFunctional - Set the given basis functional for this dual space 114 115 Not Collective 116 117 Input Parameters: 118 + sp - the `PetscDualSpace` 119 . func - the basis index 120 - q - the basis functional 121 122 Level: intermediate 123 124 Note: 125 The quadrature will be reweighted so that it has unit volume. 126 127 .seealso: `PETSCDUALSPACESIMPLE`, `PetscDualSpace`, `PetscDualSpaceSimpleSetDimension()` 128 @*/ 129 PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q) 130 { 131 PetscFunctionBegin; 132 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 133 PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace, PetscInt, PetscQuadrature), (sp, func, q)); 134 PetscFunctionReturn(PETSC_SUCCESS); 135 } 136 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 } 154 155 /*MC 156 PETSCDUALSPACESIMPLE = "simple" - A `PetscDualSpaceType` that encapsulates a dual space of functionals provided with `PetscDualSpaceSimpleSetFunctional()` 157 158 Level: intermediate 159 160 Developer Note: 161 It is not clear this has a good name 162 163 .seealso: `PetscDualSpace`, `PetscDualSpaceSimpleSetFunctional()`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()` 164 M*/ 165 166 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp) 167 { 168 PetscDualSpace_Simple *s; 169 170 PetscFunctionBegin; 171 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 172 PetscCall(PetscNew(&s)); 173 sp->data = s; 174 175 s->dim = 0; 176 s->numDof = NULL; 177 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 } 183