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