1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 3 static PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscOptionItems *PetscOptionsObject,PetscSpace sp) 4 { 5 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 6 PetscErrorCode ierr; 7 8 PetscFunctionBegin; 9 ierr = PetscOptionsHead(PetscOptionsObject,"PetscSpace polynomial options");CHKERRQ(ierr); 10 ierr = PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);CHKERRQ(ierr); 11 ierr = PetscOptionsTail();CHKERRQ(ierr); 12 PetscFunctionReturn(0); 13 } 14 15 static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v) 16 { 17 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 18 PetscErrorCode ierr; 19 20 PetscFunctionBegin; 21 ierr = PetscViewerASCIIPrintf(v, "%s space of degree %D\n", poly->tensor ? "Tensor polynomial" : "Polynomial", sp->degree);CHKERRQ(ierr); 22 PetscFunctionReturn(0); 23 } 24 25 static PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer) 26 { 27 PetscBool iascii; 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 32 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 33 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 34 if (iascii) {ierr = PetscSpacePolynomialView_Ascii(sp, viewer);CHKERRQ(ierr);} 35 PetscFunctionReturn(0); 36 } 37 38 static PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp) 39 { 40 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 41 PetscErrorCode ierr; 42 43 PetscFunctionBegin; 44 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", NULL);CHKERRQ(ierr); 45 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", NULL);CHKERRQ(ierr); 46 if (poly->subspaces) { 47 PetscInt d; 48 49 for (d = 0; d < sp->Nv; ++d) { 50 ierr = PetscSpaceDestroy(&poly->subspaces[d]);CHKERRQ(ierr); 51 } 52 } 53 ierr = PetscFree(poly->subspaces);CHKERRQ(ierr); 54 ierr = PetscFree(poly);CHKERRQ(ierr); 55 PetscFunctionReturn(0); 56 } 57 58 static PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp) 59 { 60 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 61 PetscErrorCode ierr; 62 63 PetscFunctionBegin; 64 if (poly->setupCalled) PetscFunctionReturn(0); 65 if (sp->Nv <= 1) { 66 poly->tensor = PETSC_FALSE; 67 } 68 if (sp->Nc != 1) { 69 PetscInt Nc = sp->Nc; 70 PetscBool tensor = poly->tensor; 71 PetscInt Nv = sp->Nv; 72 PetscInt degree = sp->degree; 73 PetscSpace subsp; 74 75 ierr = PetscSpaceSetType(sp, PETSCSPACESUM);CHKERRQ(ierr); 76 ierr = PetscSpaceSumSetNumSubspaces(sp, Nc);CHKERRQ(ierr); 77 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp);CHKERRQ(ierr); 78 ierr = PetscSpaceSetType(subsp, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 79 ierr = PetscSpaceSetDegree(subsp, degree, PETSC_DETERMINE);CHKERRQ(ierr); 80 ierr = PetscSpaceSetNumComponents(subsp, 1);CHKERRQ(ierr); 81 ierr = PetscSpaceSetNumVariables(subsp, Nv);CHKERRQ(ierr); 82 ierr = PetscSpacePolynomialSetTensor(subsp, tensor);CHKERRQ(ierr); 83 ierr = PetscSpaceSetUp(subsp);CHKERRQ(ierr); 84 for (PetscInt i = 0; i < Nc; i++) { 85 ierr = PetscSpaceSumSetSubspace(sp, i, subsp);CHKERRQ(ierr); 86 } 87 ierr = PetscSpaceDestroy(&subsp);CHKERRQ(ierr); 88 ierr = PetscSpaceSetUp(sp);CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 if (poly->tensor) { 92 sp->maxDegree = PETSC_DETERMINE; 93 ierr = PetscSpaceSetType(sp, PETSCSPACETENSOR);CHKERRQ(ierr); 94 ierr = PetscSpaceSetUp(sp);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 if (sp->degree < 0) SETERRQ1(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %D invalid\n", sp->degree); 98 sp->maxDegree = sp->degree; 99 poly->setupCalled = PETSC_TRUE; 100 PetscFunctionReturn(0); 101 } 102 103 static PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim) 104 { 105 PetscInt deg = sp->degree; 106 PetscInt n = sp->Nv; 107 PetscErrorCode ierr; 108 109 PetscFunctionBegin; 110 ierr = PetscDTBinomialInt(n + deg, n, dim);CHKERRQ(ierr); 111 *dim *= sp->Nc; 112 PetscFunctionReturn(0); 113 } 114 115 static PetscErrorCode CoordinateBasis(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt jet, PetscInt Njet, PetscReal pScalar[]) 116 { 117 PetscErrorCode ierr; 118 119 PetscFunctionBegin; 120 ierr = PetscArrayzero(pScalar, (1 + dim) * Njet * npoints);CHKERRQ(ierr); 121 for (PetscInt b = 0; b < 1 + dim; b++) { 122 for (PetscInt j = 0; j < PetscMin(1 + dim, Njet); j++) { 123 if (j == 0) { 124 if (b == 0) { 125 for (PetscInt pt = 0; pt < npoints; pt++) { 126 pScalar[b * Njet * npoints + j * npoints + pt] = 1.; 127 } 128 } else { 129 for (PetscInt pt = 0; pt < npoints; pt++) { 130 pScalar[b * Njet * npoints + j * npoints + pt] = points[pt * dim + (b-1)]; 131 } 132 } 133 } else if (j == b) { 134 for (PetscInt pt = 0; pt < npoints; pt++) { 135 pScalar[b * Njet * npoints + j * npoints + pt] = 1.; 136 } 137 } 138 } 139 } 140 PetscFunctionReturn(0); 141 } 142 143 static PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[]) 144 { 145 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 146 DM dm = sp->dm; 147 PetscInt dim = sp->Nv; 148 PetscInt Nb, jet, Njet; 149 PetscReal *pScalar; 150 PetscErrorCode ierr; 151 152 PetscFunctionBegin; 153 if (!poly->setupCalled) { 154 ierr = PetscSpaceSetUp(sp);CHKERRQ(ierr); 155 ierr = PetscSpaceEvaluate(sp, npoints, points, B, D, H);CHKERRQ(ierr); 156 PetscFunctionReturn(0); 157 } 158 if (poly->tensor || sp->Nc != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "tensor and multicomponent spaces should have been converted"); 159 ierr = PetscDTBinomialInt(dim + sp->degree, dim, &Nb);CHKERRQ(ierr); 160 if (H) { 161 jet = 2; 162 } else if (D) { 163 jet = 1; 164 } else { 165 jet = 0; 166 } 167 ierr = PetscDTBinomialInt(dim + jet, dim, &Njet);CHKERRQ(ierr); 168 ierr = DMGetWorkArray(dm, Nb * Njet * npoints, MPIU_REAL, &pScalar);CHKERRQ(ierr); 169 // Why are we handling the case degree == 1 specially? Because we don't want numerical noise when we evaluate hat 170 // functions at the vertices of a simplex, which happens when we invert the Vandermonde matrix of the PKD basis. 171 // We don't make any promise about which basis is used. 172 if (sp->degree == 1) { 173 ierr = CoordinateBasis(dim, npoints, points, jet, Njet, pScalar);CHKERRQ(ierr); 174 } else { 175 ierr = PetscDTPKDEvalJet(dim, npoints, points, sp->degree, jet, pScalar);CHKERRQ(ierr); 176 } 177 if (B) { 178 PetscInt p_strl = Nb; 179 PetscInt b_strl = 1; 180 181 PetscInt b_strr = Njet*npoints; 182 PetscInt p_strr = 1; 183 184 ierr = PetscArrayzero(B, npoints * Nb);CHKERRQ(ierr); 185 for (PetscInt b = 0; b < Nb; b++) { 186 for (PetscInt p = 0; p < npoints; p++) { 187 B[p*p_strl + b*b_strl] = pScalar[b*b_strr + p*p_strr]; 188 } 189 } 190 } 191 if (D) { 192 PetscInt p_strl = dim*Nb; 193 PetscInt b_strl = dim; 194 PetscInt d_strl = 1; 195 196 PetscInt b_strr = Njet*npoints; 197 PetscInt d_strr = npoints; 198 PetscInt p_strr = 1; 199 200 ierr = PetscArrayzero(D, npoints * Nb * dim);CHKERRQ(ierr); 201 for (PetscInt d = 0; d < dim; d++) { 202 for (PetscInt b = 0; b < Nb; b++) { 203 for (PetscInt p = 0; p < npoints; p++) { 204 D[p*p_strl + b*b_strl + d*d_strl] = pScalar[b*b_strr + (1+d)*d_strr + p*p_strr]; 205 } 206 } 207 } 208 } 209 if (H) { 210 PetscInt p_strl = dim*dim*Nb; 211 PetscInt b_strl = dim*dim; 212 PetscInt d1_strl = dim; 213 PetscInt d2_strl = 1; 214 215 PetscInt b_strr = Njet*npoints; 216 PetscInt j_strr = npoints; 217 PetscInt p_strr = 1; 218 219 PetscInt *derivs; 220 ierr = PetscCalloc1(dim, &derivs);CHKERRQ(ierr); 221 ierr = PetscArrayzero(H, npoints * Nb * dim * dim);CHKERRQ(ierr); 222 for (PetscInt d1 = 0; d1 < dim; d1++) { 223 for (PetscInt d2 = 0; d2 < dim; d2++) { 224 PetscInt j; 225 derivs[d1]++; 226 derivs[d2]++; 227 ierr = PetscDTGradedOrderToIndex(dim, derivs, &j);CHKERRQ(ierr); 228 derivs[d1]--; 229 derivs[d2]--; 230 for (PetscInt b = 0; b < Nb; b++) { 231 for (PetscInt p = 0; p < npoints; p++) { 232 H[p*p_strl + b*b_strl + d1*d1_strl + d2*d2_strl] = pScalar[b*b_strr + j*j_strr + p*p_strr]; 233 } 234 } 235 } 236 } 237 ierr = PetscFree(derivs);CHKERRQ(ierr); 238 } 239 ierr = DMRestoreWorkArray(dm, Nb * Njet * npoints, MPIU_REAL, &pScalar);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 } 242 243 /*@ 244 PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials (the space is spanned 245 by polynomials whose degree in each variable is bounded by the given order), as opposed to polynomials (the space is 246 spanned by polynomials whose total degree---summing over all variables---is bounded by the given order). 247 248 Input Parameters: 249 + sp - the function space object 250 - tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space 251 252 Options Database: 253 . -petscspace_poly_tensor <bool> - Whether to use tensor product polynomials in higher dimension 254 255 Level: intermediate 256 257 .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables() 258 @*/ 259 PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor) 260 { 261 PetscErrorCode ierr; 262 263 PetscFunctionBegin; 264 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 265 ierr = PetscTryMethod(sp,"PetscSpacePolynomialSetTensor_C",(PetscSpace,PetscBool),(sp,tensor));CHKERRQ(ierr); 266 PetscFunctionReturn(0); 267 } 268 269 /*@ 270 PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor polynomials (the space is spanned 271 by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is 272 spanned by polynomials whose total degree---summing over all variables---is bounded by the given order). 273 274 Input Parameters: 275 . sp - the function space object 276 277 Output Parameters: 278 . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space 279 280 Level: intermediate 281 282 .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables() 283 @*/ 284 PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor) 285 { 286 PetscErrorCode ierr; 287 288 PetscFunctionBegin; 289 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 290 PetscValidPointer(tensor, 2); 291 ierr = PetscTryMethod(sp,"PetscSpacePolynomialGetTensor_C",(PetscSpace,PetscBool*),(sp,tensor));CHKERRQ(ierr); 292 PetscFunctionReturn(0); 293 } 294 295 static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor) 296 { 297 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 298 299 PetscFunctionBegin; 300 poly->tensor = tensor; 301 PetscFunctionReturn(0); 302 } 303 304 static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor) 305 { 306 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 307 308 PetscFunctionBegin; 309 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 310 PetscValidPointer(tensor, 2); 311 *tensor = poly->tensor; 312 PetscFunctionReturn(0); 313 } 314 315 static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp) 316 { 317 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 318 PetscInt Nc, dim, order; 319 PetscBool tensor; 320 PetscErrorCode ierr; 321 322 PetscFunctionBegin; 323 ierr = PetscSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr); 324 ierr = PetscSpaceGetNumVariables(sp, &dim);CHKERRQ(ierr); 325 ierr = PetscSpaceGetDegree(sp, &order, NULL);CHKERRQ(ierr); 326 ierr = PetscSpacePolynomialGetTensor(sp, &tensor);CHKERRQ(ierr); 327 if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim); 328 if (!poly->subspaces) {ierr = PetscCalloc1(dim, &poly->subspaces);CHKERRQ(ierr);} 329 if (height <= dim) { 330 if (!poly->subspaces[height-1]) { 331 PetscSpace sub; 332 const char *name; 333 334 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);CHKERRQ(ierr); 335 ierr = PetscObjectGetName((PetscObject) sp, &name);CHKERRQ(ierr); 336 ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr); 337 ierr = PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 338 ierr = PetscSpaceSetNumComponents(sub, Nc);CHKERRQ(ierr); 339 ierr = PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);CHKERRQ(ierr); 340 ierr = PetscSpaceSetNumVariables(sub, dim-height);CHKERRQ(ierr); 341 ierr = PetscSpacePolynomialSetTensor(sub, tensor);CHKERRQ(ierr); 342 ierr = PetscSpaceSetUp(sub);CHKERRQ(ierr); 343 poly->subspaces[height-1] = sub; 344 } 345 *subsp = poly->subspaces[height-1]; 346 } else { 347 *subsp = NULL; 348 } 349 PetscFunctionReturn(0); 350 } 351 352 static PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp) 353 { 354 PetscErrorCode ierr; 355 356 PetscFunctionBegin; 357 sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial; 358 sp->ops->setup = PetscSpaceSetUp_Polynomial; 359 sp->ops->view = PetscSpaceView_Polynomial; 360 sp->ops->destroy = PetscSpaceDestroy_Polynomial; 361 sp->ops->getdimension = PetscSpaceGetDimension_Polynomial; 362 sp->ops->evaluate = PetscSpaceEvaluate_Polynomial; 363 sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial; 364 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial);CHKERRQ(ierr); 365 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial);CHKERRQ(ierr); 366 PetscFunctionReturn(0); 367 } 368 369 /*MC 370 PETSCSPACEPOLYNOMIAL = "poly" - A PetscSpace object that encapsulates a polynomial space, e.g. P1 is the space of 371 linear polynomials. The space is replicated for each component. 372 373 Level: intermediate 374 375 .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType() 376 M*/ 377 378 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp) 379 { 380 PetscSpace_Poly *poly; 381 PetscErrorCode ierr; 382 383 PetscFunctionBegin; 384 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 385 ierr = PetscNewLog(sp,&poly);CHKERRQ(ierr); 386 sp->data = poly; 387 388 poly->tensor = PETSC_FALSE; 389 poly->subspaces = NULL; 390 391 ierr = PetscSpaceInitialize_Polynomial(sp);CHKERRQ(ierr); 392 PetscFunctionReturn(0); 393 } 394 395