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