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