1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petscdmplex.h> 3 4 static PetscErrorCode PetscDualSpaceSetUp_Simple(PetscDualSpace sp) { 5 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data; 6 DM dm = sp->dm; 7 PetscInt dim, pStart, pEnd; 8 PetscSection section; 9 10 PetscFunctionBegin; 11 PetscCall(DMGetDimension(dm, &dim)); 12 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 13 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, §ion)); 14 PetscCall(PetscSectionSetChart(section, pStart, pEnd)); 15 PetscCall(PetscSectionSetDof(section, pStart, s->dim)); 16 PetscCall(PetscSectionSetUp(section)); 17 sp->pointSection = section; 18 PetscFunctionReturn(0); 19 } 20 21 static PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp) { 22 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data; 23 24 PetscFunctionBegin; 25 PetscCall(PetscFree(s->numDof)); 26 PetscCall(PetscFree(s)); 27 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetDimension_C", NULL)); 28 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetFunctional_C", NULL)); 29 PetscFunctionReturn(0); 30 } 31 32 static PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace spNew) { 33 PetscInt dim, d; 34 35 PetscFunctionBegin; 36 PetscCall(PetscDualSpaceGetDimension(sp, &dim)); 37 PetscCall(PetscDualSpaceSimpleSetDimension(spNew, dim)); 38 for (d = 0; d < dim; ++d) { 39 PetscQuadrature q; 40 41 PetscCall(PetscDualSpaceGetFunctional(sp, d, &q)); 42 PetscCall(PetscDualSpaceSimpleSetFunctional(spNew, d, q)); 43 } 44 PetscFunctionReturn(0); 45 } 46 47 static PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscDualSpace sp, PetscOptionItems *PetscOptionsObject) { 48 PetscFunctionBegin; 49 PetscFunctionReturn(0); 50 } 51 52 static PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim) { 53 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *)sp->data; 54 DM dm; 55 PetscInt spatialDim, f; 56 57 PetscFunctionBegin; 58 for (f = 0; f < s->dim; ++f) PetscCall(PetscQuadratureDestroy(&sp->functional[f])); 59 PetscCall(PetscFree(sp->functional)); 60 s->dim = dim; 61 PetscCall(PetscCalloc1(s->dim, &sp->functional)); 62 PetscCall(PetscFree(s->numDof)); 63 PetscCall(PetscDualSpaceGetDM(sp, &dm)); 64 PetscCall(DMGetCoordinateDim(dm, &spatialDim)); 65 PetscCall(PetscCalloc1(spatialDim + 1, &s->numDof)); 66 s->numDof[spatialDim] = dim; 67 PetscFunctionReturn(0); 68 } 69 70 static PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q) { 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(0); 87 } 88 89 /*@ 90 PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis 91 92 Logically Collective on sp 93 94 Input Parameters: 95 + sp - the PetscDualSpace 96 - dim - the basis dimension 97 98 Level: intermediate 99 100 .seealso: `PetscDualSpaceSimpleSetFunctional()` 101 @*/ 102 PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim) { 103 PetscFunctionBegin; 104 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 105 PetscValidLogicalCollectiveInt(sp, dim, 2); 106 PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change dimension after dual space is set up"); 107 PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace, PetscInt), (sp, dim)); 108 PetscFunctionReturn(0); 109 } 110 111 /*@ 112 PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space 113 114 Not Collective 115 116 Input Parameters: 117 + sp - the PetscDualSpace 118 . f - the basis index 119 - q - the basis functional 120 121 Level: intermediate 122 123 Note: The quadrature will be reweighted so that it has unit volume. 124 125 .seealso: `PetscDualSpaceSimpleSetDimension()` 126 @*/ 127 PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q) { 128 PetscFunctionBegin; 129 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 130 PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace, PetscInt, PetscQuadrature), (sp, func, q)); 131 PetscFunctionReturn(0); 132 } 133 134 static PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp) { 135 PetscFunctionBegin; 136 sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Simple; 137 sp->ops->setup = PetscDualSpaceSetUp_Simple; 138 sp->ops->view = NULL; 139 sp->ops->destroy = PetscDualSpaceDestroy_Simple; 140 sp->ops->duplicate = PetscDualSpaceDuplicate_Simple; 141 sp->ops->createheightsubspace = NULL; 142 sp->ops->createpointsubspace = NULL; 143 sp->ops->getsymmetries = NULL; 144 sp->ops->apply = PetscDualSpaceApplyDefault; 145 sp->ops->applyall = PetscDualSpaceApplyAllDefault; 146 sp->ops->applyint = PetscDualSpaceApplyInteriorDefault; 147 sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault; 148 sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault; 149 PetscFunctionReturn(0); 150 } 151 152 /*MC 153 PETSCDUALSPACESIMPLE = "simple" - A PetscDualSpace object that encapsulates a dual space of arbitrary functionals 154 155 Level: intermediate 156 157 .seealso: `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()` 158 M*/ 159 160 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp) { 161 PetscDualSpace_Simple *s; 162 163 PetscFunctionBegin; 164 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 165 PetscCall(PetscNewLog(sp, &s)); 166 sp->data = s; 167 168 s->dim = 0; 169 s->numDof = NULL; 170 171 PetscCall(PetscDualSpaceInitialize_Simple(sp)); 172 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple)); 173 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple)); 174 PetscFunctionReturn(0); 175 } 176