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 Options Database: 415 . -petscspace_poly_tensor <bool> - Whether to use tensor product polynomials in higher dimension 416 417 Level: beginner 418 419 .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables() 420 @*/ 421 PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor) 422 { 423 PetscErrorCode ierr; 424 425 PetscFunctionBegin; 426 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 427 ierr = PetscTryMethod(sp,"PetscSpacePolynomialSetTensor_C",(PetscSpace,PetscBool),(sp,tensor));CHKERRQ(ierr); 428 PetscFunctionReturn(0); 429 } 430 431 /*@ 432 PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor polynomials (the space is spanned 433 by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is 434 spanned by polynomials whose total degree---summing over all variables---is bounded by the given order). 435 436 Input Parameters: 437 . sp - the function space object 438 439 Output Parameters: 440 . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space 441 442 Level: beginner 443 444 .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables() 445 @*/ 446 PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor) 447 { 448 PetscErrorCode ierr; 449 450 PetscFunctionBegin; 451 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 452 PetscValidPointer(tensor, 2); 453 ierr = PetscTryMethod(sp,"PetscSpacePolynomialGetTensor_C",(PetscSpace,PetscBool*),(sp,tensor));CHKERRQ(ierr); 454 PetscFunctionReturn(0); 455 } 456 457 static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor) 458 { 459 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 460 461 PetscFunctionBegin; 462 poly->tensor = tensor; 463 PetscFunctionReturn(0); 464 } 465 466 static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor) 467 { 468 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 469 470 PetscFunctionBegin; 471 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 472 PetscValidPointer(tensor, 2); 473 *tensor = poly->tensor; 474 PetscFunctionReturn(0); 475 } 476 477 static PetscErrorCode PetscSpaceGetHeightSubspace_Polynomial(PetscSpace sp, PetscInt height, PetscSpace *subsp) 478 { 479 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 480 PetscInt Nc, dim, order; 481 PetscBool tensor; 482 PetscErrorCode ierr; 483 484 PetscFunctionBegin; 485 ierr = PetscSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr); 486 ierr = PetscSpaceGetNumVariables(sp, &dim);CHKERRQ(ierr); 487 ierr = PetscSpaceGetDegree(sp, &order, NULL);CHKERRQ(ierr); 488 ierr = PetscSpacePolynomialGetTensor(sp, &tensor);CHKERRQ(ierr); 489 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);} 490 if (!poly->subspaces) {ierr = PetscCalloc1(dim, &poly->subspaces);CHKERRQ(ierr);} 491 if (height <= dim) { 492 if (!poly->subspaces[height-1]) { 493 PetscSpace sub; 494 const char *name; 495 496 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);CHKERRQ(ierr); 497 ierr = PetscObjectGetName((PetscObject) sp, &name);CHKERRQ(ierr); 498 ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr); 499 ierr = PetscSpaceSetType(sub, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 500 ierr = PetscSpaceSetNumComponents(sub, Nc);CHKERRQ(ierr); 501 ierr = PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);CHKERRQ(ierr); 502 ierr = PetscSpaceSetNumVariables(sub, dim-height);CHKERRQ(ierr); 503 ierr = PetscSpacePolynomialSetTensor(sub, tensor);CHKERRQ(ierr); 504 ierr = PetscSpaceSetUp(sub);CHKERRQ(ierr); 505 poly->subspaces[height-1] = sub; 506 } 507 *subsp = poly->subspaces[height-1]; 508 } else { 509 *subsp = NULL; 510 } 511 PetscFunctionReturn(0); 512 } 513 514 PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp) 515 { 516 PetscErrorCode ierr; 517 518 PetscFunctionBegin; 519 sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial; 520 sp->ops->setup = PetscSpaceSetUp_Polynomial; 521 sp->ops->view = PetscSpaceView_Polynomial; 522 sp->ops->destroy = PetscSpaceDestroy_Polynomial; 523 sp->ops->getdimension = PetscSpaceGetDimension_Polynomial; 524 sp->ops->evaluate = PetscSpaceEvaluate_Polynomial; 525 sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Polynomial; 526 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial);CHKERRQ(ierr); 527 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial);CHKERRQ(ierr); 528 PetscFunctionReturn(0); 529 } 530 531 /*MC 532 PETSCSPACEPOLYNOMIAL = "poly" - A PetscSpace object that encapsulates a polynomial space, e.g. P1 is the space of 533 linear polynomials. The space is replicated for each component. 534 535 Level: intermediate 536 537 .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType() 538 M*/ 539 540 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp) 541 { 542 PetscSpace_Poly *poly; 543 PetscErrorCode ierr; 544 545 PetscFunctionBegin; 546 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 547 ierr = PetscNewLog(sp,&poly);CHKERRQ(ierr); 548 sp->data = poly; 549 550 poly->symmetric = PETSC_FALSE; 551 poly->tensor = PETSC_FALSE; 552 poly->degrees = NULL; 553 poly->subspaces = NULL; 554 555 ierr = PetscSpaceInitialize_Polynomial(sp);CHKERRQ(ierr); 556 PetscFunctionReturn(0); 557 } 558 559 PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym) 560 { 561 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 562 563 PetscFunctionBegin; 564 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 565 poly->symmetric = sym; 566 PetscFunctionReturn(0); 567 } 568 569 PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym) 570 { 571 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 572 573 PetscFunctionBegin; 574 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 575 PetscValidPointer(sym, 2); 576 *sym = poly->symmetric; 577 PetscFunctionReturn(0); 578 } 579 580