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