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 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 14 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 15 ierr = PetscSectionCreate(PETSC_COMM_SELF, §ion);CHKERRQ(ierr); 16 ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr); 17 ierr = PetscSectionSetDof(section, pStart, s->dim);CHKERRQ(ierr); 18 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 19 sp->pointSection = section; 20 PetscFunctionReturn(0); 21 } 22 23 static PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp) 24 { 25 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data; 26 PetscErrorCode ierr; 27 28 PetscFunctionBegin; 29 ierr = PetscFree(s->numDof);CHKERRQ(ierr); 30 ierr = PetscFree(s);CHKERRQ(ierr); 31 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", NULL);CHKERRQ(ierr); 32 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", NULL);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 static PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace spNew) 37 { 38 PetscInt dim, d; 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 ierr = PetscDualSpaceGetDimension(sp, &dim);CHKERRQ(ierr); 43 ierr = PetscDualSpaceSimpleSetDimension(spNew, dim);CHKERRQ(ierr); 44 for (d = 0; d < dim; ++d) { 45 PetscQuadrature q; 46 47 ierr = PetscDualSpaceGetFunctional(sp, d, &q);CHKERRQ(ierr); 48 ierr = PetscDualSpaceSimpleSetFunctional(spNew, d, q);CHKERRQ(ierr); 49 } 50 PetscFunctionReturn(0); 51 } 52 53 static PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp) 54 { 55 PetscFunctionBegin; 56 PetscFunctionReturn(0); 57 } 58 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; 65 66 PetscFunctionBegin; 67 for (f = 0; f < s->dim; ++f) {ierr = PetscQuadratureDestroy(&sp->functional[f]);CHKERRQ(ierr);} 68 ierr = PetscFree(sp->functional);CHKERRQ(ierr); 69 s->dim = dim; 70 ierr = PetscCalloc1(s->dim, &sp->functional);CHKERRQ(ierr); 71 ierr = PetscFree(s->numDof);CHKERRQ(ierr); 72 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 73 ierr = DMGetCoordinateDim(dm, &spatialDim);CHKERRQ(ierr); 74 ierr = PetscCalloc1(spatialDim+1, &s->numDof);CHKERRQ(ierr); 75 s->numDof[spatialDim] = dim; 76 PetscFunctionReturn(0); 77 } 78 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; 85 86 PetscFunctionBegin; 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 ierr = PetscQuadratureDuplicate(q, &sp->functional[f]);CHKERRQ(ierr); 89 /* Reweight so that it has unit volume: Do we want to do this for Nc > 1? */ 90 ierr = PetscQuadratureGetData(sp->functional[f], NULL, &Nc, &Nq, NULL, (const PetscReal **) &weights);CHKERRQ(ierr); 91 for (c = 0; c < Nc; ++c) { 92 PetscReal vol = 0.0; 93 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 PetscFunctionReturn(0); 98 } 99 100 /*@ 101 PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis 102 103 Logically Collective on sp 104 105 Input Parameters: 106 + sp - the PetscDualSpace 107 - dim - the basis dimension 108 109 Level: intermediate 110 111 .seealso: PetscDualSpaceSimpleSetFunctional() 112 @*/ 113 PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim) 114 { 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 119 PetscValidLogicalCollectiveInt(sp, dim, 2); 120 if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change dimension after dual space is set up"); 121 ierr = PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));CHKERRQ(ierr); 122 PetscFunctionReturn(0); 123 } 124 125 /*@ 126 PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space 127 128 Not Collective 129 130 Input Parameters: 131 + sp - the PetscDualSpace 132 . f - the basis index 133 - q - the basis functional 134 135 Level: intermediate 136 137 Note: The quadrature will be reweighted so that it has unit volume. 138 139 .seealso: PetscDualSpaceSimpleSetDimension() 140 @*/ 141 PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q) 142 { 143 PetscErrorCode ierr; 144 145 PetscFunctionBegin; 146 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 147 ierr = PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 static PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp) 152 { 153 PetscFunctionBegin; 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 PetscFunctionReturn(0); 168 } 169 170 /*MC 171 PETSCDUALSPACESIMPLE = "simple" - A PetscDualSpace object that encapsulates a dual space of arbitrary functionals 172 173 Level: intermediate 174 175 .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType() 176 M*/ 177 178 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp) 179 { 180 PetscDualSpace_Simple *s; 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 185 ierr = PetscNewLog(sp,&s);CHKERRQ(ierr); 186 sp->data = s; 187 188 s->dim = 0; 189 s->numDof = NULL; 190 191 ierr = PetscDualSpaceInitialize_Simple(sp);CHKERRQ(ierr); 192 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);CHKERRQ(ierr); 193 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196