1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petscdmplex.h> 3 4 PetscErrorCode PetscDualSpaceSetUp_Simple(PetscDualSpace sp) 5 { 6 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data; 7 DM dm = sp->dm; 8 PetscInt dim; 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 13 ierr = PetscCalloc1(dim+1, &s->numDof);CHKERRQ(ierr); 14 PetscFunctionReturn(0); 15 } 16 17 PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp) 18 { 19 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data; 20 PetscErrorCode ierr; 21 22 PetscFunctionBegin; 23 ierr = PetscFree(s->numDof);CHKERRQ(ierr); 24 ierr = PetscFree(s);CHKERRQ(ierr); 25 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", NULL);CHKERRQ(ierr); 26 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", NULL);CHKERRQ(ierr); 27 PetscFunctionReturn(0); 28 } 29 30 PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace *spNew) 31 { 32 PetscInt dim, d, Nc; 33 PetscErrorCode ierr; 34 35 PetscFunctionBegin; 36 ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);CHKERRQ(ierr); 37 ierr = PetscDualSpaceSetType(*spNew, PETSCDUALSPACESIMPLE);CHKERRQ(ierr); 38 ierr = PetscDualSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr); 39 ierr = PetscDualSpaceSetNumComponents(sp, Nc);CHKERRQ(ierr); 40 ierr = PetscDualSpaceGetDimension(sp, &dim);CHKERRQ(ierr); 41 ierr = PetscDualSpaceSimpleSetDimension(*spNew, dim);CHKERRQ(ierr); 42 for (d = 0; d < dim; ++d) { 43 PetscQuadrature q; 44 45 ierr = PetscDualSpaceGetFunctional(sp, d, &q);CHKERRQ(ierr); 46 ierr = PetscDualSpaceSimpleSetFunctional(*spNew, d, q);CHKERRQ(ierr); 47 } 48 PetscFunctionReturn(0); 49 } 50 51 PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp) 52 { 53 PetscFunctionBegin; 54 PetscFunctionReturn(0); 55 } 56 57 PetscErrorCode PetscDualSpaceGetDimension_Simple(PetscDualSpace sp, PetscInt *dim) 58 { 59 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data; 60 61 PetscFunctionBegin; 62 *dim = s->dim; 63 PetscFunctionReturn(0); 64 } 65 66 PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim) 67 { 68 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data; 69 DM dm; 70 PetscInt spatialDim, f; 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 for (f = 0; f < s->dim; ++f) {ierr = PetscQuadratureDestroy(&sp->functional[f]);CHKERRQ(ierr);} 75 ierr = PetscFree(sp->functional);CHKERRQ(ierr); 76 s->dim = dim; 77 ierr = PetscCalloc1(s->dim, &sp->functional);CHKERRQ(ierr); 78 ierr = PetscFree(s->numDof);CHKERRQ(ierr); 79 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 80 ierr = DMGetCoordinateDim(dm, &spatialDim);CHKERRQ(ierr); 81 ierr = PetscCalloc1(spatialDim+1, &s->numDof);CHKERRQ(ierr); 82 s->numDof[spatialDim] = dim; 83 PetscFunctionReturn(0); 84 } 85 86 PetscErrorCode PetscDualSpaceGetNumDof_Simple(PetscDualSpace sp, const PetscInt **numDof) 87 { 88 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data; 89 90 PetscFunctionBegin; 91 *numDof = s->numDof; 92 PetscFunctionReturn(0); 93 } 94 95 PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q) 96 { 97 PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data; 98 PetscReal *weights; 99 PetscInt Nc, c, Nq, p; 100 PetscErrorCode ierr; 101 102 PetscFunctionBegin; 103 if ((f < 0) || (f >= s->dim)) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %d not in [0, %d)", f, s->dim); 104 ierr = PetscQuadratureDuplicate(q, &sp->functional[f]);CHKERRQ(ierr); 105 /* Reweight so that it has unit volume: Do we want to do this for Nc > 1? */ 106 ierr = PetscQuadratureGetData(sp->functional[f], NULL, &Nc, &Nq, NULL, (const PetscReal **) &weights);CHKERRQ(ierr); 107 for (c = 0; c < Nc; ++c) { 108 PetscReal vol = 0.0; 109 110 for (p = 0; p < Nq; ++p) vol += weights[p*Nc+c]; 111 for (p = 0; p < Nq; ++p) weights[p*Nc+c] /= (vol == 0.0 ? 1.0 : vol); 112 } 113 PetscFunctionReturn(0); 114 } 115 116 /*@ 117 PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis 118 119 Logically Collective on PetscDualSpace 120 121 Input Parameters: 122 + sp - the PetscDualSpace 123 - dim - the basis dimension 124 125 Level: intermediate 126 127 .keywords: PetscDualSpace, dimension 128 .seealso: PetscDualSpaceSimpleSetFunctional() 129 @*/ 130 PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim) 131 { 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 136 PetscValidLogicalCollectiveInt(sp, dim, 2); 137 ierr = PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));CHKERRQ(ierr); 138 PetscFunctionReturn(0); 139 } 140 141 /*@ 142 PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space 143 144 Not Collective 145 146 Input Parameters: 147 + sp - the PetscDualSpace 148 . f - the basis index 149 - q - the basis functional 150 151 Level: intermediate 152 153 Note: The quadrature will be reweighted so that it has unit volume. 154 155 .keywords: PetscDualSpace, functional 156 .seealso: PetscDualSpaceSimpleSetDimension() 157 @*/ 158 PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q) 159 { 160 PetscErrorCode ierr; 161 162 PetscFunctionBegin; 163 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 164 ierr = PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));CHKERRQ(ierr); 165 PetscFunctionReturn(0); 166 } 167 168 PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp) 169 { 170 PetscFunctionBegin; 171 sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Simple; 172 sp->ops->setup = PetscDualSpaceSetUp_Simple; 173 sp->ops->view = NULL; 174 sp->ops->destroy = PetscDualSpaceDestroy_Simple; 175 sp->ops->duplicate = PetscDualSpaceDuplicate_Simple; 176 sp->ops->getdimension = PetscDualSpaceGetDimension_Simple; 177 sp->ops->getnumdof = PetscDualSpaceGetNumDof_Simple; 178 sp->ops->getheightsubspace = NULL; 179 sp->ops->getsymmetries = NULL; 180 sp->ops->apply = PetscDualSpaceApplyDefault; 181 sp->ops->applyall = PetscDualSpaceApplyAllDefault; 182 sp->ops->createallpoints = PetscDualSpaceCreateAllPointsDefault; 183 PetscFunctionReturn(0); 184 } 185 186 /*MC 187 PETSCDUALSPACESIMPLE = "simple" - A PetscDualSpace object that encapsulates a dual space of arbitrary functionals 188 189 Level: intermediate 190 191 .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType() 192 M*/ 193 194 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp) 195 { 196 PetscDualSpace_Simple *s; 197 PetscErrorCode ierr; 198 199 PetscFunctionBegin; 200 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 201 ierr = PetscNewLog(sp,&s);CHKERRQ(ierr); 202 sp->data = s; 203 204 s->dim = 0; 205 s->numDof = NULL; 206 207 ierr = PetscDualSpaceInitialize_Simple(sp);CHKERRQ(ierr); 208 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);CHKERRQ(ierr); 209 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);CHKERRQ(ierr); 210 PetscFunctionReturn(0); 211 } 212 213 214 215