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