1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 3 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_sym", "Use only symmetric polynomials", "PetscSpacePolynomialSetSymmetric", poly->symmetric, &poly->symmetric, NULL);CHKERRQ(ierr); 11 ierr = PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);CHKERRQ(ierr); 12 ierr = PetscOptionsTail();CHKERRQ(ierr); 13 PetscFunctionReturn(0); 14 } 15 16 static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v) 17 { 18 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 19 PetscErrorCode ierr; 20 21 PetscFunctionBegin; 22 ierr = PetscViewerASCIIPrintf(v, "%s space of degree %D\n", poly->tensor ? "Tensor polynomial" : "Polynomial", sp->degree);CHKERRQ(ierr); 23 PetscFunctionReturn(0); 24 } 25 26 PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer) 27 { 28 PetscBool iascii; 29 PetscErrorCode ierr; 30 31 PetscFunctionBegin; 32 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 33 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 34 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 35 if (iascii) {ierr = PetscSpacePolynomialView_Ascii(sp, viewer);CHKERRQ(ierr);} 36 PetscFunctionReturn(0); 37 } 38 39 PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp) 40 { 41 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 42 PetscInt ndegree = sp->degree+1; 43 PetscInt deg; 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 if (poly->setupCalled) PetscFunctionReturn(0); 48 ierr = PetscMalloc1(ndegree, &poly->degrees);CHKERRQ(ierr); 49 for (deg = 0; deg < ndegree; ++deg) poly->degrees[deg] = deg; 50 if (poly->tensor) { 51 sp->maxDegree = sp->degree + PetscMax(sp->Nv - 1,0); 52 } else { 53 sp->maxDegree = sp->degree; 54 } 55 poly->setupCalled = PETSC_TRUE; 56 PetscFunctionReturn(0); 57 } 58 59 PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp) 60 { 61 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 62 PetscErrorCode ierr; 63 64 PetscFunctionBegin; 65 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", NULL);CHKERRQ(ierr); 66 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", NULL);CHKERRQ(ierr); 67 ierr = PetscFree(poly->degrees);CHKERRQ(ierr); 68 if (poly->subspaces) { 69 PetscInt d; 70 71 for (d = 0; d < sp->Nv; ++d) { 72 ierr = PetscSpaceDestroy(&poly->subspaces[d]);CHKERRQ(ierr); 73 } 74 } 75 ierr = PetscFree(poly->subspaces);CHKERRQ(ierr); 76 ierr = PetscFree(poly);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 /* We treat the space as a tensor product of scalar polynomial spaces, so the dimension is multiplied by Nc */ 81 PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim) 82 { 83 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 84 PetscInt deg = sp->degree; 85 PetscInt n = sp->Nv, i; 86 PetscReal D = 1.0; 87 88 PetscFunctionBegin; 89 if (poly->tensor) { 90 *dim = 1; 91 for (i = 0; i < n; ++i) *dim *= (deg+1); 92 } else { 93 for (i = 1; i <= n; ++i) { 94 D *= ((PetscReal) (deg+i))/i; 95 } 96 *dim = (PetscInt) (D + 0.5); 97 } 98 *dim *= sp->Nc; 99 PetscFunctionReturn(0); 100 } 101 102 /* 103 LatticePoint_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to 'sum'. 104 105 Input Parameters: 106 + len - The length of the tuple 107 . sum - The sum of all entries in the tuple 108 - ind - The current multi-index of the tuple, initialized to the 0 tuple 109 110 Output Parameter: 111 + ind - The multi-index of the tuple, -1 indicates the iteration has terminated 112 . tup - A tuple of len integers addig to sum 113 114 Level: developer 115 116 .seealso: 117 */ 118 static PetscErrorCode LatticePoint_Internal(PetscInt len, PetscInt sum, PetscInt ind[], PetscInt tup[]) 119 { 120 PetscInt i; 121 PetscErrorCode ierr; 122 123 PetscFunctionBegin; 124 if (len == 1) { 125 ind[0] = -1; 126 tup[0] = sum; 127 } else if (sum == 0) { 128 for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;} 129 } else { 130 tup[0] = sum - ind[0]; 131 ierr = LatticePoint_Internal(len-1, ind[0], &ind[1], &tup[1]);CHKERRQ(ierr); 132 if (ind[1] < 0) { 133 if (ind[0] == sum) {ind[0] = -1;} 134 else {ind[1] = 0; ++ind[0];} 135 } 136 } 137 PetscFunctionReturn(0); 138 } 139 140 /* 141 TensorPoint_Internal - Returns all tuples of size 'len' with nonnegative integers that are less than 'max'. 142 143 Input Parameters: 144 + len - The length of the tuple 145 . max - The max for all entries in the tuple 146 - ind - The current multi-index of the tuple, initialized to the 0 tuple 147 148 Output Parameter: 149 + ind - The multi-index of the tuple, -1 indicates the iteration has terminated 150 . tup - A tuple of len integers less than max 151 152 Level: developer 153 154 .seealso: 155 */ 156 static PetscErrorCode TensorPoint_Internal(PetscInt len, PetscInt max, PetscInt ind[], PetscInt tup[]) 157 { 158 PetscInt i; 159 PetscErrorCode ierr; 160 161 PetscFunctionBegin; 162 if (len == 1) { 163 tup[0] = ind[0]++; 164 ind[0] = ind[0] >= max ? -1 : ind[0]; 165 } else if (max == 0) { 166 for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;} 167 } else { 168 tup[0] = ind[0]; 169 ierr = TensorPoint_Internal(len-1, max, &ind[1], &tup[1]);CHKERRQ(ierr); 170 if (ind[1] < 0) { 171 ind[1] = 0; 172 if (ind[0] == max-1) {ind[0] = -1;} 173 else {++ind[0];} 174 } 175 } 176 PetscFunctionReturn(0); 177 } 178 179 /* 180 p in [0, npoints), i in [0, pdim), c in [0, Nc) 181 182 B[p][i][c] = B[p][i_scalar][c][c] 183 */ 184 PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[]) 185 { 186 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 187 DM dm = sp->dm; 188 PetscInt Nc = sp->Nc; 189 PetscInt ndegree = sp->degree+1; 190 PetscInt *degrees = poly->degrees; 191 PetscInt dim = sp->Nv; 192 PetscReal *lpoints, *tmp, *LB, *LD, *LH; 193 PetscInt *ind, *tup; 194 PetscInt c, pdim, d, e, der, der2, i, p, deg, o; 195 PetscErrorCode ierr; 196 197 PetscFunctionBegin; 198 ierr = PetscSpaceGetDimension(sp, &pdim);CHKERRQ(ierr); 199 pdim /= Nc; 200 ierr = DMGetWorkArray(dm, npoints, MPIU_REAL, &lpoints);CHKERRQ(ierr); 201 ierr = DMGetWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);CHKERRQ(ierr); 202 if (B || D || H) {ierr = DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);CHKERRQ(ierr);} 203 if (D || H) {ierr = DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);CHKERRQ(ierr);} 204 if (H) {ierr = DMGetWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);CHKERRQ(ierr);} 205 for (d = 0; d < dim; ++d) { 206 for (p = 0; p < npoints; ++p) { 207 lpoints[p] = points[p*dim+d]; 208 } 209 ierr = PetscDTLegendreEval(npoints, lpoints, ndegree, degrees, tmp, &tmp[1*npoints*ndegree], &tmp[2*npoints*ndegree]);CHKERRQ(ierr); 210 /* LB, LD, LH (ndegree * dim x npoints) */ 211 for (deg = 0; deg < ndegree; ++deg) { 212 for (p = 0; p < npoints; ++p) { 213 if (B || D || H) LB[(deg*dim + d)*npoints + p] = tmp[(0*npoints + p)*ndegree+deg]; 214 if (D || H) LD[(deg*dim + d)*npoints + p] = tmp[(1*npoints + p)*ndegree+deg]; 215 if (H) LH[(deg*dim + d)*npoints + p] = tmp[(2*npoints + p)*ndegree+deg]; 216 } 217 } 218 } 219 /* Multiply by A (pdim x ndegree * dim) */ 220 ierr = PetscMalloc2(dim,&ind,dim,&tup);CHKERRQ(ierr); 221 if (B) { 222 /* B (npoints x pdim x Nc) */ 223 ierr = PetscMemzero(B, npoints*pdim*Nc*Nc * sizeof(PetscReal));CHKERRQ(ierr); 224 if (poly->tensor) { 225 i = 0; 226 ierr = PetscMemzero(ind, dim * sizeof(PetscInt));CHKERRQ(ierr); 227 while (ind[0] >= 0) { 228 ierr = TensorPoint_Internal(dim, sp->degree+1, ind, tup);CHKERRQ(ierr); 229 for (p = 0; p < npoints; ++p) { 230 B[(p*pdim + i)*Nc*Nc] = 1.0; 231 for (d = 0; d < dim; ++d) { 232 B[(p*pdim + i)*Nc*Nc] *= LB[(tup[d]*dim + d)*npoints + p]; 233 } 234 } 235 ++i; 236 } 237 } else { 238 i = 0; 239 for (o = 0; o <= sp->degree; ++o) { 240 ierr = PetscMemzero(ind, dim * sizeof(PetscInt));CHKERRQ(ierr); 241 while (ind[0] >= 0) { 242 ierr = LatticePoint_Internal(dim, o, ind, tup);CHKERRQ(ierr); 243 for (p = 0; p < npoints; ++p) { 244 B[(p*pdim + i)*Nc*Nc] = 1.0; 245 for (d = 0; d < dim; ++d) { 246 B[(p*pdim + i)*Nc*Nc] *= LB[(tup[d]*dim + d)*npoints + p]; 247 } 248 } 249 ++i; 250 } 251 } 252 } 253 /* Make direct sum basis for multicomponent space */ 254 for (p = 0; p < npoints; ++p) { 255 for (i = 0; i < pdim; ++i) { 256 for (c = 1; c < Nc; ++c) { 257 B[(p*pdim*Nc + i*Nc + c)*Nc + c] = B[(p*pdim + i)*Nc*Nc]; 258 } 259 } 260 } 261 } 262 if (D) { 263 /* D (npoints x pdim x Nc x dim) */ 264 ierr = PetscMemzero(D, npoints*pdim*Nc*Nc*dim * sizeof(PetscReal));CHKERRQ(ierr); 265 if (poly->tensor) { 266 i = 0; 267 ierr = PetscMemzero(ind, dim * sizeof(PetscInt));CHKERRQ(ierr); 268 while (ind[0] >= 0) { 269 ierr = TensorPoint_Internal(dim, sp->degree+1, ind, tup);CHKERRQ(ierr); 270 for (p = 0; p < npoints; ++p) { 271 for (der = 0; der < dim; ++der) { 272 D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0; 273 for (d = 0; d < dim; ++d) { 274 if (d == der) { 275 D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p]; 276 } else { 277 D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p]; 278 } 279 } 280 } 281 } 282 ++i; 283 } 284 } else { 285 i = 0; 286 for (o = 0; o <= sp->degree; ++o) { 287 ierr = PetscMemzero(ind, dim * sizeof(PetscInt));CHKERRQ(ierr); 288 while (ind[0] >= 0) { 289 ierr = LatticePoint_Internal(dim, o, ind, tup);CHKERRQ(ierr); 290 for (p = 0; p < npoints; ++p) { 291 for (der = 0; der < dim; ++der) { 292 D[(p*pdim + i)*Nc*Nc*dim + der] = 1.0; 293 for (d = 0; d < dim; ++d) { 294 if (d == der) { 295 D[(p*pdim + i)*Nc*Nc*dim + der] *= LD[(tup[d]*dim + d)*npoints + p]; 296 } else { 297 D[(p*pdim + i)*Nc*Nc*dim + der] *= LB[(tup[d]*dim + d)*npoints + p]; 298 } 299 } 300 } 301 } 302 ++i; 303 } 304 } 305 } 306 /* Make direct sum basis for multicomponent space */ 307 for (p = 0; p < npoints; ++p) { 308 for (i = 0; i < pdim; ++i) { 309 for (c = 1; c < Nc; ++c) { 310 for (d = 0; d < dim; ++d) { 311 D[((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d] = D[(p*pdim + i)*Nc*Nc*dim + d]; 312 } 313 } 314 } 315 } 316 } 317 if (H) { 318 /* H (npoints x pdim x Nc x Nc x dim x dim) */ 319 ierr = PetscMemzero(H, npoints*pdim*Nc*Nc*dim*dim * sizeof(PetscReal));CHKERRQ(ierr); 320 if (poly->tensor) { 321 i = 0; 322 ierr = PetscMemzero(ind, dim * sizeof(PetscInt));CHKERRQ(ierr); 323 while (ind[0] >= 0) { 324 ierr = TensorPoint_Internal(dim, sp->degree+1, ind, tup);CHKERRQ(ierr); 325 for (p = 0; p < npoints; ++p) { 326 for (der = 0; der < dim; ++der) { 327 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] = 1.0; 328 for (d = 0; d < dim; ++d) { 329 if (d == der) { 330 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LH[(tup[d]*dim + d)*npoints + p]; 331 } else { 332 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der] *= LB[(tup[d]*dim + d)*npoints + p]; 333 } 334 } 335 for (der2 = der + 1; der2 < dim; ++der2) { 336 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0; 337 for (d = 0; d < dim; ++d) { 338 if (d == der || d == der2) { 339 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p]; 340 } else { 341 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p]; 342 } 343 } 344 H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2]; 345 } 346 } 347 } 348 ++i; 349 } 350 } else { 351 i = 0; 352 for (o = 0; o <= sp->degree; ++o) { 353 ierr = PetscMemzero(ind, dim * sizeof(PetscInt));CHKERRQ(ierr); 354 while (ind[0] >= 0) { 355 ierr = LatticePoint_Internal(dim, o, ind, tup);CHKERRQ(ierr); 356 for (p = 0; p < npoints; ++p) { 357 for (der = 0; der < dim; ++der) { 358 H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] = 1.0; 359 for (d = 0; d < dim; ++d) { 360 if (d == der) { 361 H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LH[(tup[d]*dim + d)*npoints + p]; 362 } else { 363 H[((p*pdim + i)*Nc*Nc*dim + der)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p]; 364 } 365 } 366 for (der2 = der + 1; der2 < dim; ++der2) { 367 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] = 1.0; 368 for (d = 0; d < dim; ++d) { 369 if (d == der || d == der2) { 370 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LD[(tup[d]*dim + d)*npoints + p]; 371 } else { 372 H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2] *= LB[(tup[d]*dim + d)*npoints + p]; 373 } 374 } 375 H[((p*pdim + i)*Nc*Nc*dim + der2) * dim + der] = H[((p*pdim + i)*Nc*Nc*dim + der) * dim + der2]; 376 } 377 } 378 } 379 ++i; 380 } 381 } 382 } 383 /* Make direct sum basis for multicomponent space */ 384 for (p = 0; p < npoints; ++p) { 385 for (i = 0; i < pdim; ++i) { 386 for (c = 1; c < Nc; ++c) { 387 for (d = 0; d < dim; ++d) { 388 for (e = 0; e < dim; ++e) { 389 H[(((p*pdim*Nc + i*Nc + c)*Nc + c)*dim + d)*dim + e] = H[((p*pdim + i)*Nc*Nc*dim + d)*dim + e]; 390 } 391 } 392 } 393 } 394 } 395 } 396 ierr = PetscFree2(ind,tup);CHKERRQ(ierr); 397 if (H) {ierr = DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LH);CHKERRQ(ierr);} 398 if (D || H) {ierr = DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LD);CHKERRQ(ierr);} 399 if (B || D || H) {ierr = DMRestoreWorkArray(dm, npoints*dim*ndegree, MPIU_REAL, &LB);CHKERRQ(ierr);} 400 ierr = DMRestoreWorkArray(dm, npoints*ndegree*3, MPIU_REAL, &tmp);CHKERRQ(ierr); 401 ierr = DMRestoreWorkArray(dm, npoints, MPIU_REAL, &lpoints);CHKERRQ(ierr); 402 PetscFunctionReturn(0); 403 } 404 405 /*@ 406 PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials (the space is spanned 407 by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is 408 spanned by polynomials whose total degree---summing over all variables---is bounded by the given order). 409 410 Input Parameters: 411 + sp - the function space object 412 - tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space 413 414 Level: beginner 415 416 .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables() 417 @*/ 418 PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor) 419 { 420 PetscErrorCode ierr; 421 422 PetscFunctionBegin; 423 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 424 ierr = PetscTryMethod(sp,"PetscSpacePolynomialSetTensor_C",(PetscSpace,PetscBool),(sp,tensor));CHKERRQ(ierr); 425 PetscFunctionReturn(0); 426 } 427 428 /*@ 429 PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor polynomials (the space is spanned 430 by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is 431 spanned by polynomials whose total degree---summing over all variables---is bounded by the given order). 432 433 Input Parameters: 434 . sp - the function space object 435 436 Output Parameters: 437 . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space 438 439 Level: beginner 440 441 .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables() 442 @*/ 443 PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor) 444 { 445 PetscErrorCode ierr; 446 447 PetscFunctionBegin; 448 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 449 PetscValidPointer(tensor, 2); 450 ierr = PetscTryMethod(sp,"PetscSpacePolynomialGetTensor_C",(PetscSpace,PetscBool*),(sp,tensor));CHKERRQ(ierr); 451 PetscFunctionReturn(0); 452 } 453 454 static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor) 455 { 456 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 457 458 PetscFunctionBegin; 459 poly->tensor = tensor; 460 PetscFunctionReturn(0); 461 } 462 463 static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor) 464 { 465 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 466 467 PetscFunctionBegin; 468 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 469 PetscValidPointer(tensor, 2); 470 *tensor = poly->tensor; 471 PetscFunctionReturn(0); 472 } 473 474 static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp) 475 { 476 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 477 PetscInt Nc, dim, order; 478 PetscBool tensor; 479 PetscErrorCode ierr; 480 481 PetscFunctionBegin; 482 ierr = PetscSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr); 483 ierr = PetscSpaceGetNumVariables(sp, &dim);CHKERRQ(ierr); 484 ierr = PetscSpaceGetDegree(sp, &order, NULL);CHKERRQ(ierr); 485 ierr = PetscSpacePolynomialGetTensor(sp, &tensor);CHKERRQ(ierr); 486 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);} 487 if (!poly->subspaces) {ierr = PetscCalloc1(dim, &poly->subspaces);CHKERRQ(ierr);} 488 if (height <= dim) { 489 if (!poly->subspaces[height-1]) { 490 PetscSpace sub; 491 const char *name; 492 493 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);CHKERRQ(ierr); 494 ierr = PetscObjectGetName((PetscObject) sp, &name);CHKERRQ(ierr); 495 ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr); 496 ierr = PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 497 ierr = PetscSpaceSetNumComponents(sub, Nc);CHKERRQ(ierr); 498 ierr = PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);CHKERRQ(ierr); 499 ierr = PetscSpaceSetNumVariables(sub, dim-height);CHKERRQ(ierr); 500 ierr = PetscSpacePolynomialSetTensor(sub, tensor);CHKERRQ(ierr); 501 ierr = PetscSpaceSetUp(sub);CHKERRQ(ierr); 502 poly->subspaces[height-1] = sub; 503 } 504 *subsp = poly->subspaces[height-1]; 505 } else { 506 *subsp = NULL; 507 } 508 PetscFunctionReturn(0); 509 } 510 511 PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp) 512 { 513 PetscErrorCode ierr; 514 515 PetscFunctionBegin; 516 sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial; 517 sp->ops->setup = PetscSpaceSetUp_Polynomial; 518 sp->ops->view = PetscSpaceView_Polynomial; 519 sp->ops->destroy = PetscSpaceDestroy_Polynomial; 520 sp->ops->getdimension = PetscSpaceGetDimension_Polynomial; 521 sp->ops->evaluate = PetscSpaceEvaluate_Polynomial; 522 sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial; 523 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial);CHKERRQ(ierr); 524 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial);CHKERRQ(ierr); 525 PetscFunctionReturn(0); 526 } 527 528 /*MC 529 PETSCSPACEPOLYNOMIAL = "poly" - A PetscSpace object that encapsulates a polynomial space, e.g. P1 is the space of 530 linear polynomials. The space is replicated for each component. 531 532 Level: intermediate 533 534 .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType() 535 M*/ 536 537 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp) 538 { 539 PetscSpace_Poly *poly; 540 PetscErrorCode ierr; 541 542 PetscFunctionBegin; 543 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 544 ierr = PetscNewLog(sp,&poly);CHKERRQ(ierr); 545 sp->data = poly; 546 547 poly->symmetric = PETSC_FALSE; 548 poly->tensor = PETSC_FALSE; 549 poly->degrees = NULL; 550 poly->subspaces = NULL; 551 552 ierr = PetscSpaceInitialize_Polynomial(sp);CHKERRQ(ierr); 553 PetscFunctionReturn(0); 554 } 555 556 PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym) 557 { 558 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 559 560 PetscFunctionBegin; 561 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 562 poly->symmetric = sym; 563 PetscFunctionReturn(0); 564 } 565 566 PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym) 567 { 568 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 569 570 PetscFunctionBegin; 571 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 572 PetscValidPointer(sym, 2); 573 *sym = poly->symmetric; 574 PetscFunctionReturn(0); 575 } 576 577