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; 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 static 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 static 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 static PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp) 52 { 53 PetscFunctionBegin; 54 PetscFunctionReturn(0); 55 } 56 57 static 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 static 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 static 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 static 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 sp 120 121 Input Parameters: 122 + sp - the PetscDualSpace 123 - dim - the basis dimension 124 125 Level: intermediate 126 127 .seealso: PetscDualSpaceSimpleSetFunctional() 128 @*/ 129 PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim) 130 { 131 PetscErrorCode ierr; 132 133 PetscFunctionBegin; 134 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 135 PetscValidLogicalCollectiveInt(sp, dim, 2); 136 ierr = PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 140 /*@ 141 PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space 142 143 Not Collective 144 145 Input Parameters: 146 + sp - the PetscDualSpace 147 . f - the basis index 148 - q - the basis functional 149 150 Level: intermediate 151 152 Note: The quadrature will be reweighted so that it has unit volume. 153 154 .seealso: PetscDualSpaceSimpleSetDimension() 155 @*/ 156 PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q) 157 { 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 162 ierr = PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));CHKERRQ(ierr); 163 PetscFunctionReturn(0); 164 } 165 166 static PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp) 167 { 168 PetscFunctionBegin; 169 sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Simple; 170 sp->ops->setup = PetscDualSpaceSetUp_Simple; 171 sp->ops->view = NULL; 172 sp->ops->destroy = PetscDualSpaceDestroy_Simple; 173 sp->ops->duplicate = PetscDualSpaceDuplicate_Simple; 174 sp->ops->getdimension = PetscDualSpaceGetDimension_Simple; 175 sp->ops->getnumdof = PetscDualSpaceGetNumDof_Simple; 176 sp->ops->getheightsubspace = NULL; 177 sp->ops->getsymmetries = NULL; 178 sp->ops->apply = PetscDualSpaceApplyDefault; 179 sp->ops->applyall = PetscDualSpaceApplyAllDefault; 180 sp->ops->createallpoints = PetscDualSpaceCreateAllPointsDefault; 181 PetscFunctionReturn(0); 182 } 183 184 /*MC 185 PETSCDUALSPACESIMPLE = "simple" - A PetscDualSpace object that encapsulates a dual space of arbitrary functionals 186 187 Level: intermediate 188 189 .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType() 190 M*/ 191 192 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp) 193 { 194 PetscDualSpace_Simple *s; 195 PetscErrorCode ierr; 196 197 PetscFunctionBegin; 198 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 199 ierr = PetscNewLog(sp,&s);CHKERRQ(ierr); 200 sp->data = s; 201 202 s->dim = 0; 203 s->numDof = NULL; 204 205 ierr = PetscDualSpaceInitialize_Simple(sp);CHKERRQ(ierr); 206 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);CHKERRQ(ierr); 207 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);CHKERRQ(ierr); 208 PetscFunctionReturn(0); 209 } 210