1 /* Discretization tools */ 2 3 #include <petscdt.h> /*I "petscdt.h" I*/ 4 #include <petscblaslapack.h> 5 #include <petsc/private/petscimpl.h> 6 #include <petsc/private/dtimpl.h> 7 #include <petsc/private/petscfeimpl.h> /* For CoordinatesRefToReal() */ 8 #include <petscviewer.h> 9 #include <petscdmplex.h> 10 #include <petscdmshell.h> 11 12 #if defined(PETSC_HAVE_MPFR) 13 #include <mpfr.h> 14 #endif 15 16 const char *const PetscDTNodeTypes_shifted[] = {"default", "gaussjacobi", "equispaced", "tanhsinh", "PETSCDTNODES_", NULL}; 17 const char *const *const PetscDTNodeTypes = PetscDTNodeTypes_shifted + 1; 18 19 const char *const PetscDTSimplexQuadratureTypes_shifted[] = {"default", "conic", "minsym", "PETSCDTSIMPLEXQUAD_", NULL}; 20 const char *const *const PetscDTSimplexQuadratureTypes = PetscDTSimplexQuadratureTypes_shifted + 1; 21 22 static PetscBool GolubWelschCite = PETSC_FALSE; 23 const char GolubWelschCitation[] = "@article{GolubWelsch1969,\n" 24 " author = {Golub and Welsch},\n" 25 " title = {Calculation of Quadrature Rules},\n" 26 " journal = {Math. Comp.},\n" 27 " volume = {23},\n" 28 " number = {106},\n" 29 " pages = {221--230},\n" 30 " year = {1969}\n}\n"; 31 32 /* Numerical tests in src/dm/dt/tests/ex1.c show that when computing the nodes and weights of Gauss-Jacobi 33 quadrature rules: 34 35 - in double precision, Newton's method and Golub & Welsch both work for moderate degrees (< 100), 36 - in single precision, Newton's method starts producing incorrect roots around n = 15, but 37 the weights from Golub & Welsch become a problem before then: they produces errors 38 in computing the Jacobi-polynomial Gram matrix around n = 6. 39 40 So we default to Newton's method (required fewer dependencies) */ 41 PetscBool PetscDTGaussQuadratureNewton_Internal = PETSC_TRUE; 42 43 PetscClassId PETSCQUADRATURE_CLASSID = 0; 44 45 /*@ 46 PetscQuadratureCreate - Create a `PetscQuadrature` object 47 48 Collective 49 50 Input Parameter: 51 . comm - The communicator for the `PetscQuadrature` object 52 53 Output Parameter: 54 . q - The `PetscQuadrature` object 55 56 Level: beginner 57 58 .seealso: `PetscQuadrature`, `Petscquadraturedestroy()`, `PetscQuadratureGetData()` 59 @*/ 60 PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q) 61 { 62 PetscFunctionBegin; 63 PetscAssertPointer(q, 2); 64 PetscCall(DMInitializePackage()); 65 PetscCall(PetscHeaderCreate(*q, PETSCQUADRATURE_CLASSID, "PetscQuadrature", "Quadrature", "DT", comm, PetscQuadratureDestroy, PetscQuadratureView)); 66 (*q)->ct = DM_POLYTOPE_UNKNOWN; 67 (*q)->dim = -1; 68 (*q)->Nc = 1; 69 (*q)->order = -1; 70 (*q)->numPoints = 0; 71 (*q)->points = NULL; 72 (*q)->weights = NULL; 73 PetscFunctionReturn(PETSC_SUCCESS); 74 } 75 76 /*@ 77 PetscQuadratureDuplicate - Create a deep copy of the `PetscQuadrature` object 78 79 Collective 80 81 Input Parameter: 82 . q - The `PetscQuadrature` object 83 84 Output Parameter: 85 . r - The new `PetscQuadrature` object 86 87 Level: beginner 88 89 .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureDestroy()`, `PetscQuadratureGetData()` 90 @*/ 91 PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r) 92 { 93 DMPolytopeType ct; 94 PetscInt order, dim, Nc, Nq; 95 const PetscReal *points, *weights; 96 PetscReal *p, *w; 97 98 PetscFunctionBegin; 99 PetscAssertPointer(q, 1); 100 PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)q), r)); 101 PetscCall(PetscQuadratureGetCellType(q, &ct)); 102 PetscCall(PetscQuadratureSetCellType(*r, ct)); 103 PetscCall(PetscQuadratureGetOrder(q, &order)); 104 PetscCall(PetscQuadratureSetOrder(*r, order)); 105 PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights)); 106 PetscCall(PetscMalloc1(Nq * dim, &p)); 107 PetscCall(PetscMalloc1(Nq * Nc, &w)); 108 PetscCall(PetscArraycpy(p, points, Nq * dim)); 109 PetscCall(PetscArraycpy(w, weights, Nc * Nq)); 110 PetscCall(PetscQuadratureSetData(*r, dim, Nc, Nq, p, w)); 111 PetscFunctionReturn(PETSC_SUCCESS); 112 } 113 114 /*@ 115 PetscQuadratureDestroy - Destroys a `PetscQuadrature` object 116 117 Collective 118 119 Input Parameter: 120 . q - The `PetscQuadrature` object 121 122 Level: beginner 123 124 .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()` 125 @*/ 126 PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q) 127 { 128 PetscFunctionBegin; 129 if (!*q) PetscFunctionReturn(PETSC_SUCCESS); 130 PetscValidHeaderSpecific(*q, PETSCQUADRATURE_CLASSID, 1); 131 if (--((PetscObject)*q)->refct > 0) { 132 *q = NULL; 133 PetscFunctionReturn(PETSC_SUCCESS); 134 } 135 PetscCall(PetscFree((*q)->points)); 136 PetscCall(PetscFree((*q)->weights)); 137 PetscCall(PetscHeaderDestroy(q)); 138 PetscFunctionReturn(PETSC_SUCCESS); 139 } 140 141 /*@ 142 PetscQuadratureGetCellType - Return the cell type of the integration domain 143 144 Not Collective 145 146 Input Parameter: 147 . q - The `PetscQuadrature` object 148 149 Output Parameter: 150 . ct - The cell type of the integration domain 151 152 Level: intermediate 153 154 .seealso: `PetscQuadrature`, `PetscQuadratureSetCellType()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()` 155 @*/ 156 PetscErrorCode PetscQuadratureGetCellType(PetscQuadrature q, DMPolytopeType *ct) 157 { 158 PetscFunctionBegin; 159 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 160 PetscAssertPointer(ct, 2); 161 *ct = q->ct; 162 PetscFunctionReturn(PETSC_SUCCESS); 163 } 164 165 /*@ 166 PetscQuadratureSetCellType - Set the cell type of the integration domain 167 168 Not Collective 169 170 Input Parameters: 171 + q - The `PetscQuadrature` object 172 - ct - The cell type of the integration domain 173 174 Level: intermediate 175 176 .seealso: `PetscQuadrature`, `PetscQuadratureGetCellType()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()` 177 @*/ 178 PetscErrorCode PetscQuadratureSetCellType(PetscQuadrature q, DMPolytopeType ct) 179 { 180 PetscFunctionBegin; 181 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 182 q->ct = ct; 183 PetscFunctionReturn(PETSC_SUCCESS); 184 } 185 186 /*@ 187 PetscQuadratureGetOrder - Return the order of the method in the `PetscQuadrature` 188 189 Not Collective 190 191 Input Parameter: 192 . q - The `PetscQuadrature` object 193 194 Output Parameter: 195 . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated 196 197 Level: intermediate 198 199 .seealso: `PetscQuadrature`, `PetscQuadratureSetOrder()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()` 200 @*/ 201 PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order) 202 { 203 PetscFunctionBegin; 204 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 205 PetscAssertPointer(order, 2); 206 *order = q->order; 207 PetscFunctionReturn(PETSC_SUCCESS); 208 } 209 210 /*@ 211 PetscQuadratureSetOrder - Set the order of the method in the `PetscQuadrature` 212 213 Not Collective 214 215 Input Parameters: 216 + q - The `PetscQuadrature` object 217 - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated 218 219 Level: intermediate 220 221 .seealso: `PetscQuadrature`, `PetscQuadratureGetOrder()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()` 222 @*/ 223 PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order) 224 { 225 PetscFunctionBegin; 226 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 227 q->order = order; 228 PetscFunctionReturn(PETSC_SUCCESS); 229 } 230 231 /*@ 232 PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated 233 234 Not Collective 235 236 Input Parameter: 237 . q - The `PetscQuadrature` object 238 239 Output Parameter: 240 . Nc - The number of components 241 242 Level: intermediate 243 244 Note: 245 We are performing an integral $\int f(x) w(x) dx$, where both $f$ and $w$ (the weight) have `Nc` components. 246 247 .seealso: `PetscQuadrature`, `PetscQuadratureSetNumComponents()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()` 248 @*/ 249 PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc) 250 { 251 PetscFunctionBegin; 252 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 253 PetscAssertPointer(Nc, 2); 254 *Nc = q->Nc; 255 PetscFunctionReturn(PETSC_SUCCESS); 256 } 257 258 /*@ 259 PetscQuadratureSetNumComponents - Sets the number of components for functions to be integrated 260 261 Not Collective 262 263 Input Parameters: 264 + q - The `PetscQuadrature` object 265 - Nc - The number of components 266 267 Level: intermediate 268 269 Note: 270 We are performing an integral $\int f(x) w(x) dx$, where both $f$ and $w$ (the weight) have `Nc` components. 271 272 .seealso: `PetscQuadrature`, `PetscQuadratureGetNumComponents()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()` 273 @*/ 274 PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc) 275 { 276 PetscFunctionBegin; 277 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 278 q->Nc = Nc; 279 PetscFunctionReturn(PETSC_SUCCESS); 280 } 281 282 /*@C 283 PetscQuadratureGetData - Returns the data defining the `PetscQuadrature` 284 285 Not Collective 286 287 Input Parameter: 288 . q - The `PetscQuadrature` object 289 290 Output Parameters: 291 + dim - The spatial dimension 292 . Nc - The number of components 293 . npoints - The number of quadrature points 294 . points - The coordinates of each quadrature point 295 - weights - The weight of each quadrature point 296 297 Level: intermediate 298 299 Note: 300 All output arguments are optional, pass `NULL` for any argument not required 301 302 Fortran Note: 303 Call `PetscQuadratureRestoreData()` when you are done with the data 304 305 .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureSetData()` 306 @*/ 307 PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PeOp PetscInt *dim, PeOp PetscInt *Nc, PeOp PetscInt *npoints, PeOp const PetscReal *points[], PeOp const PetscReal *weights[]) 308 { 309 PetscFunctionBegin; 310 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 311 if (dim) { 312 PetscAssertPointer(dim, 2); 313 *dim = q->dim; 314 } 315 if (Nc) { 316 PetscAssertPointer(Nc, 3); 317 *Nc = q->Nc; 318 } 319 if (npoints) { 320 PetscAssertPointer(npoints, 4); 321 *npoints = q->numPoints; 322 } 323 if (points) { 324 PetscAssertPointer(points, 5); 325 *points = q->points; 326 } 327 if (weights) { 328 PetscAssertPointer(weights, 6); 329 *weights = q->weights; 330 } 331 PetscFunctionReturn(PETSC_SUCCESS); 332 } 333 334 /*@ 335 PetscQuadratureEqual - determine whether two quadratures are equivalent 336 337 Input Parameters: 338 + A - A `PetscQuadrature` object 339 - B - Another `PetscQuadrature` object 340 341 Output Parameter: 342 . equal - `PETSC_TRUE` if the quadratures are the same 343 344 Level: intermediate 345 346 .seealso: `PetscQuadrature`, `PetscQuadratureCreate()` 347 @*/ 348 PetscErrorCode PetscQuadratureEqual(PetscQuadrature A, PetscQuadrature B, PetscBool *equal) 349 { 350 PetscFunctionBegin; 351 PetscValidHeaderSpecific(A, PETSCQUADRATURE_CLASSID, 1); 352 PetscValidHeaderSpecific(B, PETSCQUADRATURE_CLASSID, 2); 353 PetscAssertPointer(equal, 3); 354 *equal = PETSC_FALSE; 355 if (A->ct != B->ct || A->dim != B->dim || A->Nc != B->Nc || A->order != B->order || A->numPoints != B->numPoints) PetscFunctionReturn(PETSC_SUCCESS); 356 for (PetscInt i = 0; i < A->numPoints * A->dim; i++) { 357 if (PetscAbsReal(A->points[i] - B->points[i]) > PETSC_SMALL) PetscFunctionReturn(PETSC_SUCCESS); 358 } 359 if (!A->weights && !B->weights) { 360 *equal = PETSC_TRUE; 361 PetscFunctionReturn(PETSC_SUCCESS); 362 } 363 if (A->weights && B->weights) { 364 for (PetscInt i = 0; i < A->numPoints; i++) { 365 if (PetscAbsReal(A->weights[i] - B->weights[i]) > PETSC_SMALL) PetscFunctionReturn(PETSC_SUCCESS); 366 } 367 *equal = PETSC_TRUE; 368 } 369 PetscFunctionReturn(PETSC_SUCCESS); 370 } 371 372 static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[]) 373 { 374 PetscScalar *Js, *Jinvs; 375 PetscInt i, j, k; 376 PetscBLASInt bm, bn, info; 377 378 PetscFunctionBegin; 379 if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS); 380 PetscCall(PetscBLASIntCast(m, &bm)); 381 PetscCall(PetscBLASIntCast(n, &bn)); 382 #if defined(PETSC_USE_COMPLEX) 383 PetscCall(PetscMalloc2(m * n, &Js, m * n, &Jinvs)); 384 for (i = 0; i < m * n; i++) Js[i] = J[i]; 385 #else 386 Js = (PetscReal *)J; 387 Jinvs = Jinv; 388 #endif 389 if (m == n) { 390 PetscBLASInt *pivots; 391 PetscScalar *W; 392 393 PetscCall(PetscMalloc2(m, &pivots, m, &W)); 394 395 PetscCall(PetscArraycpy(Jinvs, Js, m * m)); 396 PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info)); 397 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscBLASInt_FMT, info); 398 PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &info)); 399 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscBLASInt_FMT, info); 400 PetscCall(PetscFree2(pivots, W)); 401 } else if (m < n) { 402 PetscScalar *JJT; 403 PetscBLASInt *pivots; 404 PetscScalar *W; 405 406 PetscCall(PetscMalloc1(m * m, &JJT)); 407 PetscCall(PetscMalloc2(m, &pivots, m, &W)); 408 for (i = 0; i < m; i++) { 409 for (j = 0; j < m; j++) { 410 PetscScalar val = 0.; 411 412 for (k = 0; k < n; k++) val += Js[i * n + k] * Js[j * n + k]; 413 JJT[i * m + j] = val; 414 } 415 } 416 417 PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info)); 418 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscBLASInt_FMT, info); 419 PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info)); 420 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscBLASInt_FMT, info); 421 for (i = 0; i < n; i++) { 422 for (j = 0; j < m; j++) { 423 PetscScalar val = 0.; 424 425 for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j]; 426 Jinvs[i * m + j] = val; 427 } 428 } 429 PetscCall(PetscFree2(pivots, W)); 430 PetscCall(PetscFree(JJT)); 431 } else { 432 PetscScalar *JTJ; 433 PetscBLASInt *pivots; 434 PetscScalar *W; 435 436 PetscCall(PetscMalloc1(n * n, &JTJ)); 437 PetscCall(PetscMalloc2(n, &pivots, n, &W)); 438 for (i = 0; i < n; i++) { 439 for (j = 0; j < n; j++) { 440 PetscScalar val = 0.; 441 442 for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j]; 443 JTJ[i * n + j] = val; 444 } 445 } 446 447 PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bn, pivots, &info)); 448 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscBLASInt_FMT, info); 449 PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info)); 450 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscBLASInt_FMT, info); 451 for (i = 0; i < n; i++) { 452 for (j = 0; j < m; j++) { 453 PetscScalar val = 0.; 454 455 for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k]; 456 Jinvs[i * m + j] = val; 457 } 458 } 459 PetscCall(PetscFree2(pivots, W)); 460 PetscCall(PetscFree(JTJ)); 461 } 462 #if defined(PETSC_USE_COMPLEX) 463 for (i = 0; i < m * n; i++) Jinv[i] = PetscRealPart(Jinvs[i]); 464 PetscCall(PetscFree2(Js, Jinvs)); 465 #endif 466 PetscFunctionReturn(PETSC_SUCCESS); 467 } 468 469 /*@ 470 PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation. 471 472 Collective 473 474 Input Parameters: 475 + q - the quadrature functional 476 . imageDim - the dimension of the image of the transformation 477 . origin - a point in the original space 478 . originImage - the image of the origin under the transformation 479 . J - the Jacobian of the image: an [imageDim x dim] matrix in row major order 480 - formDegree - transform the quadrature weights as k-forms of this form degree (if the number of components is a multiple of (dim choose `formDegree`), 481 it is assumed that they represent multiple k-forms) [see `PetscDTAltVPullback()` for interpretation of `formDegree`] 482 483 Output Parameter: 484 . Jinvstarq - a quadrature rule where each point is the image of a point in the original quadrature rule, and where the k-form weights have 485 been pulled-back by the pseudoinverse of `J` to the k-form weights in the image space. 486 487 Level: intermediate 488 489 Note: 490 The new quadrature rule will have a different number of components if spaces have different dimensions. 491 For example, pushing a 2-form forward from a two dimensional space to a three dimensional space changes the number of components from 1 to 3. 492 493 .seealso: `PetscQuadrature`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()` 494 @*/ 495 PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq) 496 { 497 PetscInt dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c; 498 const PetscReal *points; 499 const PetscReal *weights; 500 PetscReal *imagePoints, *imageWeights; 501 PetscReal *Jinv; 502 PetscReal *Jinvstar; 503 504 PetscFunctionBegin; 505 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 506 PetscCheck(imageDim >= PetscAbsInt(formDegree), PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Cannot represent a %" PetscInt_FMT "-form in %" PetscInt_FMT " dimensions", PetscAbsInt(formDegree), imageDim); 507 PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights)); 508 PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize)); 509 PetscCheck(Nc % formSize == 0, PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Number of components %" PetscInt_FMT " is not a multiple of formSize %" PetscInt_FMT, Nc, formSize); 510 Ncopies = Nc / formSize; 511 PetscCall(PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize)); 512 imageNc = Ncopies * imageFormSize; 513 PetscCall(PetscMalloc1(Npoints * imageDim, &imagePoints)); 514 PetscCall(PetscMalloc1(Npoints * imageNc, &imageWeights)); 515 PetscCall(PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar)); 516 PetscCall(PetscDTJacobianInverse_Internal(imageDim, dim, J, Jinv)); 517 PetscCall(PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar)); 518 for (pt = 0; pt < Npoints; pt++) { 519 const PetscReal *point = PetscSafePointerPlusOffset(points, pt * dim); 520 PetscReal *imagePoint = &imagePoints[pt * imageDim]; 521 522 for (i = 0; i < imageDim; i++) { 523 PetscReal val = originImage[i]; 524 525 for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]); 526 imagePoint[i] = val; 527 } 528 for (c = 0; c < Ncopies; c++) { 529 const PetscReal *form = &weights[pt * Nc + c * formSize]; 530 PetscReal *imageForm = &imageWeights[pt * imageNc + c * imageFormSize]; 531 532 for (i = 0; i < imageFormSize; i++) { 533 PetscReal val = 0.; 534 535 for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j]; 536 imageForm[i] = val; 537 } 538 } 539 } 540 PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq)); 541 PetscCall(PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights)); 542 PetscCall(PetscFree2(Jinv, Jinvstar)); 543 PetscFunctionReturn(PETSC_SUCCESS); 544 } 545 546 /*@C 547 PetscQuadratureSetData - Sets the data defining the quadrature 548 549 Not Collective 550 551 Input Parameters: 552 + q - The `PetscQuadrature` object 553 . dim - The spatial dimension 554 . Nc - The number of components 555 . npoints - The number of quadrature points 556 . points - The coordinates of each quadrature point 557 - weights - The weight of each quadrature point 558 559 Level: intermediate 560 561 Note: 562 `q` owns the references to points and weights, so they must be allocated using `PetscMalloc()` and the user should not free them. 563 564 .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()` 565 @*/ 566 PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[]) PeNSS 567 { 568 PetscFunctionBegin; 569 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 570 if (dim >= 0) q->dim = dim; 571 if (Nc >= 0) q->Nc = Nc; 572 if (npoints >= 0) q->numPoints = npoints; 573 if (points) { 574 PetscAssertPointer(points, 5); 575 q->points = points; 576 } 577 if (weights) { 578 PetscAssertPointer(weights, 6); 579 q->weights = weights; 580 } 581 PetscFunctionReturn(PETSC_SUCCESS); 582 } 583 584 static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v) 585 { 586 PetscInt q, d, c; 587 PetscViewerFormat format; 588 589 PetscFunctionBegin; 590 if (quad->Nc > 1) 591 PetscCall(PetscViewerASCIIPrintf(v, "Quadrature on a %s of order %" PetscInt_FMT " on %" PetscInt_FMT " points (dim %" PetscInt_FMT ") with %" PetscInt_FMT " components\n", DMPolytopeTypes[quad->ct], quad->order, quad->numPoints, quad->dim, quad->Nc)); 592 else PetscCall(PetscViewerASCIIPrintf(v, "Quadrature on a %s of order %" PetscInt_FMT " on %" PetscInt_FMT " points (dim %" PetscInt_FMT ")\n", DMPolytopeTypes[quad->ct], quad->order, quad->numPoints, quad->dim)); 593 PetscCall(PetscViewerGetFormat(v, &format)); 594 if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS); 595 for (q = 0; q < quad->numPoints; ++q) { 596 PetscCall(PetscViewerASCIIPrintf(v, "p%" PetscInt_FMT " (", q)); 597 PetscCall(PetscViewerASCIIUseTabs(v, PETSC_FALSE)); 598 for (d = 0; d < quad->dim; ++d) { 599 if (d) PetscCall(PetscViewerASCIIPrintf(v, ", ")); 600 PetscCall(PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q * quad->dim + d])); 601 } 602 PetscCall(PetscViewerASCIIPrintf(v, ") ")); 603 if (quad->Nc > 1) PetscCall(PetscViewerASCIIPrintf(v, "w%" PetscInt_FMT " (", q)); 604 for (c = 0; c < quad->Nc; ++c) { 605 if (c) PetscCall(PetscViewerASCIIPrintf(v, ", ")); 606 PetscCall(PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q * quad->Nc + c])); 607 } 608 if (quad->Nc > 1) PetscCall(PetscViewerASCIIPrintf(v, ")")); 609 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 610 PetscCall(PetscViewerASCIIUseTabs(v, PETSC_TRUE)); 611 } 612 PetscFunctionReturn(PETSC_SUCCESS); 613 } 614 615 /*@ 616 PetscQuadratureView - View a `PetscQuadrature` object 617 618 Collective 619 620 Input Parameters: 621 + quad - The `PetscQuadrature` object 622 - viewer - The `PetscViewer` object 623 624 Level: beginner 625 626 .seealso: `PetscQuadrature`, `PetscViewer`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()` 627 @*/ 628 PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer) 629 { 630 PetscBool iascii; 631 632 PetscFunctionBegin; 633 PetscValidHeader(quad, 1); 634 if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 635 if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)quad), &viewer)); 636 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 637 PetscCall(PetscViewerASCIIPushTab(viewer)); 638 if (iascii) PetscCall(PetscQuadratureView_Ascii(quad, viewer)); 639 PetscCall(PetscViewerASCIIPopTab(viewer)); 640 PetscFunctionReturn(PETSC_SUCCESS); 641 } 642 643 /*@C 644 PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement 645 646 Not Collective; No Fortran Support 647 648 Input Parameters: 649 + q - The original `PetscQuadrature` 650 . numSubelements - The number of subelements the original element is divided into 651 . v0 - An array of the initial points for each subelement 652 - jac - An array of the Jacobian mappings from the reference to each subelement 653 654 Output Parameter: 655 . qref - The dimension 656 657 Level: intermediate 658 659 Note: 660 Together `v0` and `jac` define an affine mapping from the original reference element to each subelement 661 662 .seealso: `PetscQuadrature`, `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()` 663 @*/ 664 PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref) 665 { 666 DMPolytopeType ct; 667 const PetscReal *points, *weights; 668 PetscReal *pointsRef, *weightsRef; 669 PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e; 670 671 PetscFunctionBegin; 672 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 673 PetscAssertPointer(v0, 3); 674 PetscAssertPointer(jac, 4); 675 PetscAssertPointer(qref, 5); 676 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, qref)); 677 PetscCall(PetscQuadratureGetCellType(q, &ct)); 678 PetscCall(PetscQuadratureGetOrder(q, &order)); 679 PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights)); 680 npointsRef = npoints * numSubelements; 681 PetscCall(PetscMalloc1(npointsRef * dim, &pointsRef)); 682 PetscCall(PetscMalloc1(npointsRef * Nc, &weightsRef)); 683 for (c = 0; c < numSubelements; ++c) { 684 for (p = 0; p < npoints; ++p) { 685 for (d = 0; d < dim; ++d) { 686 pointsRef[(c * npoints + p) * dim + d] = v0[c * dim + d]; 687 for (e = 0; e < dim; ++e) pointsRef[(c * npoints + p) * dim + d] += jac[(c * dim + d) * dim + e] * (points[p * dim + e] + 1.0); 688 } 689 /* Could also use detJ here */ 690 for (cp = 0; cp < Nc; ++cp) weightsRef[(c * npoints + p) * Nc + cp] = weights[p * Nc + cp] / numSubelements; 691 } 692 } 693 PetscCall(PetscQuadratureSetCellType(*qref, ct)); 694 PetscCall(PetscQuadratureSetOrder(*qref, order)); 695 PetscCall(PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef)); 696 PetscFunctionReturn(PETSC_SUCCESS); 697 } 698 699 /* Compute the coefficients for the Jacobi polynomial recurrence, 700 701 J^{a,b}_n(x) = (cnm1 + cnm1x * x) * J^{a,b}_{n-1}(x) - cnm2 * J^{a,b}_{n-2}(x). 702 */ 703 #define PetscDTJacobiRecurrence_Internal(n, a, b, cnm1, cnm1x, cnm2) \ 704 do { \ 705 PetscReal _a = (a); \ 706 PetscReal _b = (b); \ 707 PetscReal _n = (n); \ 708 if (n == 1) { \ 709 (cnm1) = (_a - _b) * 0.5; \ 710 (cnm1x) = (_a + _b + 2.) * 0.5; \ 711 (cnm2) = 0.; \ 712 } else { \ 713 PetscReal _2n = _n + _n; \ 714 PetscReal _d = (_2n * (_n + _a + _b) * (_2n + _a + _b - 2)); \ 715 PetscReal _n1 = (_2n + _a + _b - 1.) * (_a * _a - _b * _b); \ 716 PetscReal _n1x = (_2n + _a + _b - 1.) * (_2n + _a + _b) * (_2n + _a + _b - 2); \ 717 PetscReal _n2 = 2. * ((_n + _a - 1.) * (_n + _b - 1.) * (_2n + _a + _b)); \ 718 (cnm1) = _n1 / _d; \ 719 (cnm1x) = _n1x / _d; \ 720 (cnm2) = _n2 / _d; \ 721 } \ 722 } while (0) 723 724 /*@ 725 PetscDTJacobiNorm - Compute the weighted L2 norm of a Jacobi polynomial. 726 727 $\| P^{\alpha,\beta}_n \|_{\alpha,\beta}^2 = \int_{-1}^1 (1 + x)^{\alpha} (1 - x)^{\beta} P^{\alpha,\beta}_n (x)^2 dx.$ 728 729 Input Parameters: 730 + alpha - the left exponent > -1 731 . beta - the right exponent > -1 732 - n - the polynomial degree 733 734 Output Parameter: 735 . norm - the weighted L2 norm 736 737 Level: beginner 738 739 .seealso: `PetscQuadrature`, `PetscDTJacobiEval()` 740 @*/ 741 PetscErrorCode PetscDTJacobiNorm(PetscReal alpha, PetscReal beta, PetscInt n, PetscReal *norm) 742 { 743 PetscReal twoab1; 744 PetscReal gr; 745 746 PetscFunctionBegin; 747 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent alpha %g <= -1. invalid", (double)alpha); 748 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent beta %g <= -1. invalid", (double)beta); 749 PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "n %" PetscInt_FMT " < 0 invalid", n); 750 twoab1 = PetscPowReal(2., alpha + beta + 1.); 751 #if defined(PETSC_HAVE_LGAMMA) 752 if (!n) { 753 gr = PetscExpReal(PetscLGamma(alpha + 1.) + PetscLGamma(beta + 1.) - PetscLGamma(alpha + beta + 2.)); 754 } else { 755 gr = PetscExpReal(PetscLGamma(n + alpha + 1.) + PetscLGamma(n + beta + 1.) - (PetscLGamma(n + 1.) + PetscLGamma(n + alpha + beta + 1.))) / (n + n + alpha + beta + 1.); 756 } 757 #else 758 { 759 PetscInt alphai = (PetscInt)alpha; 760 PetscInt betai = (PetscInt)beta; 761 PetscInt i; 762 763 gr = n ? (1. / (n + n + alpha + beta + 1.)) : 1.; 764 if ((PetscReal)alphai == alpha) { 765 if (!n) { 766 for (i = 0; i < alphai; i++) gr *= (i + 1.) / (beta + i + 1.); 767 gr /= (alpha + beta + 1.); 768 } else { 769 for (i = 0; i < alphai; i++) gr *= (n + i + 1.) / (n + beta + i + 1.); 770 } 771 } else if ((PetscReal)betai == beta) { 772 if (!n) { 773 for (i = 0; i < betai; i++) gr *= (i + 1.) / (alpha + i + 2.); 774 gr /= (alpha + beta + 1.); 775 } else { 776 for (i = 0; i < betai; i++) gr *= (n + i + 1.) / (n + alpha + i + 1.); 777 } 778 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable."); 779 } 780 #endif 781 *norm = PetscSqrtReal(twoab1 * gr); 782 PetscFunctionReturn(PETSC_SUCCESS); 783 } 784 785 static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p) 786 { 787 PetscReal ak, bk; 788 PetscReal abk1; 789 PetscInt i, l, maxdegree; 790 791 PetscFunctionBegin; 792 maxdegree = degrees[ndegree - 1] - k; 793 ak = a + k; 794 bk = b + k; 795 abk1 = a + b + k + 1.; 796 if (maxdegree < 0) { 797 for (i = 0; i < npoints; i++) 798 for (l = 0; l < ndegree; l++) p[i * ndegree + l] = 0.; 799 PetscFunctionReturn(PETSC_SUCCESS); 800 } 801 for (i = 0; i < npoints; i++) { 802 PetscReal pm1, pm2, x; 803 PetscReal cnm1, cnm1x, cnm2; 804 PetscInt j, m; 805 806 x = points[i]; 807 pm2 = 1.; 808 PetscDTJacobiRecurrence_Internal(1, ak, bk, cnm1, cnm1x, cnm2); 809 pm1 = (cnm1 + cnm1x * x); 810 l = 0; 811 while (l < ndegree && degrees[l] - k < 0) p[l++] = 0.; 812 while (l < ndegree && degrees[l] - k == 0) { 813 p[l] = pm2; 814 for (m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5; 815 l++; 816 } 817 while (l < ndegree && degrees[l] - k == 1) { 818 p[l] = pm1; 819 for (m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5; 820 l++; 821 } 822 for (j = 2; j <= maxdegree; j++) { 823 PetscReal pp; 824 825 PetscDTJacobiRecurrence_Internal(j, ak, bk, cnm1, cnm1x, cnm2); 826 pp = (cnm1 + cnm1x * x) * pm1 - cnm2 * pm2; 827 pm2 = pm1; 828 pm1 = pp; 829 while (l < ndegree && degrees[l] - k == j) { 830 p[l] = pp; 831 for (m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5; 832 l++; 833 } 834 } 835 p += ndegree; 836 } 837 PetscFunctionReturn(PETSC_SUCCESS); 838 } 839 840 /*@ 841 PetscDTJacobiEvalJet - Evaluate the jet (function and derivatives) of the Jacobi polynomials basis up to a given degree. 842 843 Input Parameters: 844 + alpha - the left exponent of the weight 845 . beta - the right exponetn of the weight 846 . npoints - the number of points to evaluate the polynomials at 847 . points - [npoints] array of point coordinates 848 . degree - the maximm degree polynomial space to evaluate, (degree + 1) will be evaluated total. 849 - k - the maximum derivative to evaluate in the jet, (k + 1) will be evaluated total. 850 851 Output Parameter: 852 . p - an array containing the evaluations of the Jacobi polynomials's jets on the points. the size is (degree + 1) x 853 (k + 1) x npoints, which also describes the order of the dimensions of this three-dimensional array: the first 854 (slowest varying) dimension is polynomial degree; the second dimension is derivative order; the third (fastest 855 varying) dimension is the index of the evaluation point. 856 857 Level: advanced 858 859 Notes: 860 The Jacobi polynomials with indices $\alpha$ and $\beta$ are orthogonal with respect to the 861 weighted inner product $\langle f, g \rangle = \int_{-1}^1 (1+x)^{\alpha} (1-x)^{\beta} f(x) 862 g(x) dx$. 863 864 .seealso: `PetscDTJacobiEval()`, `PetscDTPKDEvalJet()` 865 @*/ 866 PetscErrorCode PetscDTJacobiEvalJet(PetscReal alpha, PetscReal beta, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[]) 867 { 868 PetscInt i, j, l; 869 PetscInt *degrees; 870 PetscReal *psingle; 871 872 PetscFunctionBegin; 873 if (degree == 0) { 874 PetscInt zero = 0; 875 876 for (i = 0; i <= k; i++) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, 1, &zero, &p[i * npoints])); 877 PetscFunctionReturn(PETSC_SUCCESS); 878 } 879 PetscCall(PetscMalloc1(degree + 1, °rees)); 880 PetscCall(PetscMalloc1((degree + 1) * npoints, &psingle)); 881 for (i = 0; i <= degree; i++) degrees[i] = i; 882 for (i = 0; i <= k; i++) { 883 PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, degree + 1, degrees, psingle)); 884 for (j = 0; j <= degree; j++) { 885 for (l = 0; l < npoints; l++) p[(j * (k + 1) + i) * npoints + l] = psingle[l * (degree + 1) + j]; 886 } 887 } 888 PetscCall(PetscFree(psingle)); 889 PetscCall(PetscFree(degrees)); 890 PetscFunctionReturn(PETSC_SUCCESS); 891 } 892 893 /*@ 894 PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$ at a set of points 895 at points 896 897 Not Collective 898 899 Input Parameters: 900 + npoints - number of spatial points to evaluate at 901 . alpha - the left exponent > -1 902 . beta - the right exponent > -1 903 . points - array of locations to evaluate at 904 . ndegree - number of basis degrees to evaluate 905 - degrees - sorted array of degrees to evaluate 906 907 Output Parameters: 908 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or `NULL`) 909 . D - row-oriented derivative evaluation matrix (or `NULL`) 910 - D2 - row-oriented second derivative evaluation matrix (or `NULL`) 911 912 Level: intermediate 913 914 .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()` 915 @*/ 916 PetscErrorCode PetscDTJacobiEval(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PeOp PetscReal B[], PeOp PetscReal D[], PeOp PetscReal D2[]) 917 { 918 PetscFunctionBegin; 919 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1."); 920 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1."); 921 if (!npoints || !ndegree) PetscFunctionReturn(PETSC_SUCCESS); 922 if (B) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B)); 923 if (D) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D)); 924 if (D2) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2)); 925 PetscFunctionReturn(PETSC_SUCCESS); 926 } 927 928 /*@ 929 PetscDTLegendreEval - evaluate Legendre polynomials at points 930 931 Not Collective 932 933 Input Parameters: 934 + npoints - number of spatial points to evaluate at 935 . points - array of locations to evaluate at 936 . ndegree - number of basis degrees to evaluate 937 - degrees - sorted array of degrees to evaluate 938 939 Output Parameters: 940 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension `npoints`*`ndegrees`, allocated by caller) (or `NULL`) 941 . D - row-oriented derivative evaluation matrix (or `NULL`) 942 - D2 - row-oriented second derivative evaluation matrix (or `NULL`) 943 944 Level: intermediate 945 946 .seealso: `PetscDTGaussQuadrature()` 947 @*/ 948 PetscErrorCode PetscDTLegendreEval(PetscInt npoints, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PeOp PetscReal B[], PeOp PetscReal D[], PeOp PetscReal D2[]) 949 { 950 PetscFunctionBegin; 951 PetscCall(PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2)); 952 PetscFunctionReturn(PETSC_SUCCESS); 953 } 954 955 /*@ 956 PetscDTIndexToGradedOrder - convert an index into a tuple of monomial degrees in a graded order (that is, if the degree sum of tuple x is less than the degree sum of tuple y, 957 then the index of x is smaller than the index of y) 958 959 Input Parameters: 960 + len - the desired length of the degree tuple 961 - index - the index to convert: should be >= 0 962 963 Output Parameter: 964 . degtup - filled with a tuple of degrees 965 966 Level: beginner 967 968 Note: 969 For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples 970 acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the 971 last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1). 972 973 .seealso: `PetscDTGradedOrderToIndex()` 974 @*/ 975 PetscErrorCode PetscDTIndexToGradedOrder(PetscInt len, PetscInt index, PetscInt degtup[]) 976 { 977 PetscInt i, total; 978 PetscInt sum; 979 980 PetscFunctionBeginHot; 981 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 982 PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative"); 983 total = 1; 984 sum = 0; 985 while (index >= total) { 986 index -= total; 987 total = (total * (len + sum)) / (sum + 1); 988 sum++; 989 } 990 for (i = 0; i < len; i++) { 991 PetscInt c; 992 993 degtup[i] = sum; 994 for (c = 0, total = 1; c < sum; c++) { 995 /* going into the loop, total is the number of way to have a tuple of sum exactly c with length len - 1 - i */ 996 if (index < total) break; 997 index -= total; 998 total = (total * (len - 1 - i + c)) / (c + 1); 999 degtup[i]--; 1000 } 1001 sum -= degtup[i]; 1002 } 1003 PetscFunctionReturn(PETSC_SUCCESS); 1004 } 1005 1006 /*@ 1007 PetscDTGradedOrderToIndex - convert a tuple into an index in a graded order, the inverse of `PetscDTIndexToGradedOrder()`. 1008 1009 Input Parameters: 1010 + len - the length of the degree tuple 1011 - degtup - tuple with this length 1012 1013 Output Parameter: 1014 . index - index in graded order: >= 0 1015 1016 Level: beginner 1017 1018 Note: 1019 For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples 1020 acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the 1021 last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1). 1022 1023 .seealso: `PetscDTIndexToGradedOrder()` 1024 @*/ 1025 PetscErrorCode PetscDTGradedOrderToIndex(PetscInt len, const PetscInt degtup[], PetscInt *index) 1026 { 1027 PetscInt i, idx, sum, total; 1028 1029 PetscFunctionBeginHot; 1030 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 1031 for (i = 0, sum = 0; i < len; i++) sum += degtup[i]; 1032 idx = 0; 1033 total = 1; 1034 for (i = 0; i < sum; i++) { 1035 idx += total; 1036 total = (total * (len + i)) / (i + 1); 1037 } 1038 for (i = 0; i < len - 1; i++) { 1039 PetscInt c; 1040 1041 total = 1; 1042 sum -= degtup[i]; 1043 for (c = 0; c < sum; c++) { 1044 idx += total; 1045 total = (total * (len - 1 - i + c)) / (c + 1); 1046 } 1047 } 1048 *index = idx; 1049 PetscFunctionReturn(PETSC_SUCCESS); 1050 } 1051 1052 static PetscBool PKDCite = PETSC_FALSE; 1053 const char PKDCitation[] = "@article{Kirby2010,\n" 1054 " title={Singularity-free evaluation of collapsed-coordinate orthogonal polynomials},\n" 1055 " author={Kirby, Robert C},\n" 1056 " journal={ACM Transactions on Mathematical Software (TOMS)},\n" 1057 " volume={37},\n" 1058 " number={1},\n" 1059 " pages={1--16},\n" 1060 " year={2010},\n" 1061 " publisher={ACM New York, NY, USA}\n}\n"; 1062 1063 /*@ 1064 PetscDTPKDEvalJet - Evaluate the jet (function and derivatives) of the Proriol-Koornwinder-Dubiner (PKD) basis for 1065 the space of polynomials up to a given degree. 1066 1067 Input Parameters: 1068 + dim - the number of variables in the multivariate polynomials 1069 . npoints - the number of points to evaluate the polynomials at 1070 . points - [npoints x dim] array of point coordinates 1071 . degree - the degree (sum of degrees on the variables in a monomial) of the polynomial space to evaluate. There are ((dim + degree) choose dim) polynomials in this space. 1072 - k - the maximum order partial derivative to evaluate in the jet. There are (dim + k choose dim) partial derivatives 1073 in the jet. Choosing k = 0 means to evaluate just the function and no derivatives 1074 1075 Output Parameter: 1076 . p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is ((dim + degree) 1077 choose dim) x ((dim + k) choose dim) x npoints, which also describes the order of the dimensions of this 1078 three-dimensional array: the first (slowest varying) dimension is basis function index; the second dimension is jet 1079 index; the third (fastest varying) dimension is the index of the evaluation point. 1080 1081 Level: advanced 1082 1083 Notes: 1084 The PKD basis is L2-orthonormal on the biunit simplex (which is used as the reference element 1085 for finite elements in PETSc), which makes it a stable basis to use for evaluating 1086 polynomials in that domain. 1087 1088 The ordering of the basis functions, and the ordering of the derivatives in the jet, both follow the graded 1089 ordering of `PetscDTIndexToGradedOrder()` and `PetscDTGradedOrderToIndex()`. For example, in 3D, the polynomial with 1090 leading monomial x^2,y^0,z^1, which has degree tuple (2,0,1), which by `PetscDTGradedOrderToIndex()` has index 12 (it is the 13th basis function in the space); 1091 the partial derivative $\partial_x \partial_z$ has order tuple (1,0,1), appears at index 6 in the jet (it is the 7th partial derivative in the jet). 1092 1093 The implementation uses Kirby's singularity-free evaluation algorithm, <https://doi.org/10.1145/1644001.1644006>. 1094 1095 .seealso: `PetscDTGradedOrderToIndex()`, `PetscDTIndexToGradedOrder()`, `PetscDTJacobiEvalJet()` 1096 @*/ 1097 PetscErrorCode PetscDTPKDEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[]) 1098 { 1099 PetscInt degidx, kidx, d, pt; 1100 PetscInt Nk, Ndeg; 1101 PetscInt *ktup, *degtup; 1102 PetscReal *scales, initscale, scaleexp; 1103 1104 PetscFunctionBegin; 1105 PetscCall(PetscCitationsRegister(PKDCitation, &PKDCite)); 1106 PetscCall(PetscDTBinomialInt(dim + k, k, &Nk)); 1107 PetscCall(PetscDTBinomialInt(degree + dim, degree, &Ndeg)); 1108 PetscCall(PetscMalloc2(dim, °tup, dim, &ktup)); 1109 PetscCall(PetscMalloc1(Ndeg, &scales)); 1110 initscale = 1.; 1111 if (dim > 1) { 1112 PetscCall(PetscDTBinomial(dim, 2, &scaleexp)); 1113 initscale = PetscPowReal(2., scaleexp * 0.5); 1114 } 1115 for (degidx = 0; degidx < Ndeg; degidx++) { 1116 PetscInt e, i; 1117 PetscInt m1idx = -1, m2idx = -1; 1118 PetscInt n; 1119 PetscInt degsum; 1120 PetscReal alpha; 1121 PetscReal cnm1, cnm1x, cnm2; 1122 PetscReal norm; 1123 1124 PetscCall(PetscDTIndexToGradedOrder(dim, degidx, degtup)); 1125 for (d = dim - 1; d >= 0; d--) 1126 if (degtup[d]) break; 1127 if (d < 0) { /* constant is 1 everywhere, all derivatives are zero */ 1128 scales[degidx] = initscale; 1129 for (e = 0; e < dim; e++) { 1130 PetscCall(PetscDTJacobiNorm(e, 0., 0, &norm)); 1131 scales[degidx] /= norm; 1132 } 1133 for (i = 0; i < npoints; i++) p[degidx * Nk * npoints + i] = 1.; 1134 for (i = 0; i < (Nk - 1) * npoints; i++) p[(degidx * Nk + 1) * npoints + i] = 0.; 1135 continue; 1136 } 1137 n = degtup[d]; 1138 degtup[d]--; 1139 PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m1idx)); 1140 if (degtup[d] > 0) { 1141 degtup[d]--; 1142 PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m2idx)); 1143 degtup[d]++; 1144 } 1145 degtup[d]++; 1146 for (e = 0, degsum = 0; e < d; e++) degsum += degtup[e]; 1147 alpha = 2 * degsum + d; 1148 PetscDTJacobiRecurrence_Internal(n, alpha, 0., cnm1, cnm1x, cnm2); 1149 1150 scales[degidx] = initscale; 1151 for (e = 0, degsum = 0; e < dim; e++) { 1152 PetscInt f; 1153 PetscReal ealpha; 1154 PetscReal enorm; 1155 1156 ealpha = 2 * degsum + e; 1157 for (f = 0; f < degsum; f++) scales[degidx] *= 2.; 1158 PetscCall(PetscDTJacobiNorm(ealpha, 0., degtup[e], &enorm)); 1159 scales[degidx] /= enorm; 1160 degsum += degtup[e]; 1161 } 1162 1163 for (pt = 0; pt < npoints; pt++) { 1164 /* compute the multipliers */ 1165 PetscReal thetanm1, thetanm1x, thetanm2; 1166 1167 thetanm1x = dim - (d + 1) + 2. * points[pt * dim + d]; 1168 for (e = d + 1; e < dim; e++) thetanm1x += points[pt * dim + e]; 1169 thetanm1x *= 0.5; 1170 thetanm1 = (2. - (dim - (d + 1))); 1171 for (e = d + 1; e < dim; e++) thetanm1 -= points[pt * dim + e]; 1172 thetanm1 *= 0.5; 1173 thetanm2 = thetanm1 * thetanm1; 1174 1175 for (kidx = 0; kidx < Nk; kidx++) { 1176 PetscInt f; 1177 1178 PetscCall(PetscDTIndexToGradedOrder(dim, kidx, ktup)); 1179 /* first sum in the same derivative terms */ 1180 p[(degidx * Nk + kidx) * npoints + pt] = (cnm1 * thetanm1 + cnm1x * thetanm1x) * p[(m1idx * Nk + kidx) * npoints + pt]; 1181 if (m2idx >= 0) p[(degidx * Nk + kidx) * npoints + pt] -= cnm2 * thetanm2 * p[(m2idx * Nk + kidx) * npoints + pt]; 1182 1183 for (f = d; f < dim; f++) { 1184 PetscInt km1idx, mplty = ktup[f]; 1185 1186 if (!mplty) continue; 1187 ktup[f]--; 1188 PetscCall(PetscDTGradedOrderToIndex(dim, ktup, &km1idx)); 1189 1190 /* the derivative of cnm1x * thetanm1x wrt x variable f is 0.5 * cnm1x if f > d otherwise it is cnm1x */ 1191 /* the derivative of cnm1 * thetanm1 wrt x variable f is 0 if f == d, otherwise it is -0.5 * cnm1 */ 1192 /* the derivative of -cnm2 * thetanm2 wrt x variable f is 0 if f == d, otherwise it is cnm2 * thetanm1 */ 1193 if (f > d) { 1194 PetscInt f2; 1195 1196 p[(degidx * Nk + kidx) * npoints + pt] += mplty * 0.5 * (cnm1x - cnm1) * p[(m1idx * Nk + km1idx) * npoints + pt]; 1197 if (m2idx >= 0) { 1198 p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm2 * thetanm1 * p[(m2idx * Nk + km1idx) * npoints + pt]; 1199 /* second derivatives of -cnm2 * thetanm2 wrt x variable f,f2 is like - 0.5 * cnm2 */ 1200 for (f2 = f; f2 < dim; f2++) { 1201 PetscInt km2idx, mplty2 = ktup[f2]; 1202 PetscInt factor; 1203 1204 if (!mplty2) continue; 1205 ktup[f2]--; 1206 PetscCall(PetscDTGradedOrderToIndex(dim, ktup, &km2idx)); 1207 1208 factor = mplty * mplty2; 1209 if (f == f2) factor /= 2; 1210 p[(degidx * Nk + kidx) * npoints + pt] -= 0.5 * factor * cnm2 * p[(m2idx * Nk + km2idx) * npoints + pt]; 1211 ktup[f2]++; 1212 } 1213 } 1214 } else { 1215 p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm1x * p[(m1idx * Nk + km1idx) * npoints + pt]; 1216 } 1217 ktup[f]++; 1218 } 1219 } 1220 } 1221 } 1222 for (degidx = 0; degidx < Ndeg; degidx++) { 1223 PetscReal scale = scales[degidx]; 1224 PetscInt i; 1225 1226 for (i = 0; i < Nk * npoints; i++) p[degidx * Nk * npoints + i] *= scale; 1227 } 1228 PetscCall(PetscFree(scales)); 1229 PetscCall(PetscFree2(degtup, ktup)); 1230 PetscFunctionReturn(PETSC_SUCCESS); 1231 } 1232 1233 /*@ 1234 PetscDTPTrimmedSize - The size of the trimmed polynomial space of k-forms with a given degree and form degree, 1235 which can be evaluated in `PetscDTPTrimmedEvalJet()`. 1236 1237 Input Parameters: 1238 + dim - the number of variables in the multivariate polynomials 1239 . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space. 1240 - formDegree - the degree of the form 1241 1242 Output Parameter: 1243 . size - The number ((`dim` + `degree`) choose (`dim` + `formDegree`)) x ((`degree` + `formDegree` - 1) choose (`formDegree`)) 1244 1245 Level: advanced 1246 1247 .seealso: `PetscDTPTrimmedEvalJet()` 1248 @*/ 1249 PetscErrorCode PetscDTPTrimmedSize(PetscInt dim, PetscInt degree, PetscInt formDegree, PetscInt *size) 1250 { 1251 PetscInt Nrk, Nbpt; // number of trimmed polynomials 1252 1253 PetscFunctionBegin; 1254 formDegree = PetscAbsInt(formDegree); 1255 PetscCall(PetscDTBinomialInt(degree + dim, degree + formDegree, &Nbpt)); 1256 PetscCall(PetscDTBinomialInt(degree + formDegree - 1, formDegree, &Nrk)); 1257 Nbpt *= Nrk; 1258 *size = Nbpt; 1259 PetscFunctionReturn(PETSC_SUCCESS); 1260 } 1261 1262 /* there was a reference implementation based on section 4.4 of Arnold, Falk & Winther (acta numerica, 2006), but it 1263 * was inferior to this implementation */ 1264 static PetscErrorCode PetscDTPTrimmedEvalJet_Internal(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[]) 1265 { 1266 PetscInt formDegreeOrig = formDegree; 1267 PetscBool formNegative = (formDegreeOrig < 0) ? PETSC_TRUE : PETSC_FALSE; 1268 1269 PetscFunctionBegin; 1270 formDegree = PetscAbsInt(formDegreeOrig); 1271 if (formDegree == 0) { 1272 PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree, jetDegree, p)); 1273 PetscFunctionReturn(PETSC_SUCCESS); 1274 } 1275 if (formDegree == dim) { 1276 PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p)); 1277 PetscFunctionReturn(PETSC_SUCCESS); 1278 } 1279 PetscInt Nbpt; 1280 PetscCall(PetscDTPTrimmedSize(dim, degree, formDegree, &Nbpt)); 1281 PetscInt Nf; 1282 PetscCall(PetscDTBinomialInt(dim, formDegree, &Nf)); 1283 PetscInt Nk; 1284 PetscCall(PetscDTBinomialInt(dim + jetDegree, dim, &Nk)); 1285 PetscCall(PetscArrayzero(p, Nbpt * Nf * Nk * npoints)); 1286 1287 PetscInt Nbpm1; // number of scalar polynomials up to degree - 1; 1288 PetscCall(PetscDTBinomialInt(dim + degree - 1, dim, &Nbpm1)); 1289 PetscReal *p_scalar; 1290 PetscCall(PetscMalloc1(Nbpm1 * Nk * npoints, &p_scalar)); 1291 PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p_scalar)); 1292 PetscInt total = 0; 1293 // First add the full polynomials up to degree - 1 into the basis: take the scalar 1294 // and copy one for each form component 1295 for (PetscInt i = 0; i < Nbpm1; i++) { 1296 const PetscReal *src = &p_scalar[i * Nk * npoints]; 1297 for (PetscInt f = 0; f < Nf; f++) { 1298 PetscReal *dest = &p[(total++ * Nf + f) * Nk * npoints]; 1299 PetscCall(PetscArraycpy(dest, src, Nk * npoints)); 1300 } 1301 } 1302 PetscInt *form_atoms; 1303 PetscCall(PetscMalloc1(formDegree + 1, &form_atoms)); 1304 // construct the interior product pattern 1305 PetscInt (*pattern)[3]; 1306 PetscInt Nf1; // number of formDegree + 1 forms 1307 PetscCall(PetscDTBinomialInt(dim, formDegree + 1, &Nf1)); 1308 PetscInt nnz = Nf1 * (formDegree + 1); 1309 PetscCall(PetscMalloc1(Nf1 * (formDegree + 1), &pattern)); 1310 PetscCall(PetscDTAltVInteriorPattern(dim, formDegree + 1, pattern)); 1311 PetscReal centroid = (1. - dim) / (dim + 1.); 1312 PetscInt *deriv; 1313 PetscCall(PetscMalloc1(dim, &deriv)); 1314 for (PetscInt d = dim; d >= formDegree + 1; d--) { 1315 PetscInt Nfd1; // number of formDegree + 1 forms in dimension d that include dx_0 1316 // (equal to the number of formDegree forms in dimension d-1) 1317 PetscCall(PetscDTBinomialInt(d - 1, formDegree, &Nfd1)); 1318 // The number of homogeneous (degree-1) scalar polynomials in d variables 1319 PetscInt Nh; 1320 PetscCall(PetscDTBinomialInt(d - 1 + degree - 1, d - 1, &Nh)); 1321 const PetscReal *h_scalar = &p_scalar[(Nbpm1 - Nh) * Nk * npoints]; 1322 for (PetscInt b = 0; b < Nh; b++) { 1323 const PetscReal *h_s = &h_scalar[b * Nk * npoints]; 1324 for (PetscInt f = 0; f < Nfd1; f++) { 1325 // construct all formDegree+1 forms that start with dx_(dim - d) /\ ... 1326 form_atoms[0] = dim - d; 1327 PetscCall(PetscDTEnumSubset(d - 1, formDegree, f, &form_atoms[1])); 1328 for (PetscInt i = 0; i < formDegree; i++) form_atoms[1 + i] += form_atoms[0] + 1; 1329 PetscInt f_ind; // index of the resulting form 1330 PetscCall(PetscDTSubsetIndex(dim, formDegree + 1, form_atoms, &f_ind)); 1331 PetscReal *p_f = &p[total++ * Nf * Nk * npoints]; 1332 for (PetscInt nz = 0; nz < nnz; nz++) { 1333 PetscInt i = pattern[nz][0]; // formDegree component 1334 PetscInt j = pattern[nz][1]; // (formDegree + 1) component 1335 PetscInt v = pattern[nz][2]; // coordinate component 1336 PetscReal scale = v < 0 ? -1. : 1.; 1337 1338 i = formNegative ? (Nf - 1 - i) : i; 1339 scale = (formNegative && (i & 1)) ? -scale : scale; 1340 v = v < 0 ? -(v + 1) : v; 1341 if (j != f_ind) continue; 1342 PetscReal *p_i = &p_f[i * Nk * npoints]; 1343 for (PetscInt jet = 0; jet < Nk; jet++) { 1344 const PetscReal *h_jet = &h_s[jet * npoints]; 1345 PetscReal *p_jet = &p_i[jet * npoints]; 1346 1347 for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * h_jet[pt] * (points[pt * dim + v] - centroid); 1348 PetscCall(PetscDTIndexToGradedOrder(dim, jet, deriv)); 1349 deriv[v]++; 1350 PetscReal mult = deriv[v]; 1351 PetscInt l; 1352 PetscCall(PetscDTGradedOrderToIndex(dim, deriv, &l)); 1353 if (l >= Nk) continue; 1354 p_jet = &p_i[l * npoints]; 1355 for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * mult * h_jet[pt]; 1356 deriv[v]--; 1357 } 1358 } 1359 } 1360 } 1361 } 1362 PetscCheck(total == Nbpt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrectly counted P trimmed polynomials"); 1363 PetscCall(PetscFree(deriv)); 1364 PetscCall(PetscFree(pattern)); 1365 PetscCall(PetscFree(form_atoms)); 1366 PetscCall(PetscFree(p_scalar)); 1367 PetscFunctionReturn(PETSC_SUCCESS); 1368 } 1369 1370 /*@ 1371 PetscDTPTrimmedEvalJet - Evaluate the jet (function and derivatives) of a basis of the trimmed polynomial k-forms up to 1372 a given degree. 1373 1374 Input Parameters: 1375 + dim - the number of variables in the multivariate polynomials 1376 . npoints - the number of points to evaluate the polynomials at 1377 . points - [npoints x dim] array of point coordinates 1378 . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space to evaluate. 1379 There are ((dim + degree) choose (dim + formDegree)) x ((degree + formDegree - 1) choose (formDegree)) polynomials in this space. 1380 (You can use `PetscDTPTrimmedSize()` to compute this size.) 1381 . formDegree - the degree of the form 1382 - jetDegree - the maximum order partial derivative to evaluate in the jet. There are ((dim + jetDegree) choose dim) partial derivatives 1383 in the jet. Choosing jetDegree = 0 means to evaluate just the function and no derivatives 1384 1385 Output Parameter: 1386 . p - an array containing the evaluations of the PKD polynomials' jets on the points. 1387 1388 Level: advanced 1389 1390 Notes: 1391 The size of `p` is `PetscDTPTrimmedSize()` x ((dim + formDegree) choose dim) x ((dim + k) 1392 choose dim) x npoints,which also describes the order of the dimensions of this 1393 four-dimensional array\: 1394 1395 the first (slowest varying) dimension is basis function index; 1396 the second dimension is component of the form; 1397 the third dimension is jet index; 1398 the fourth (fastest varying) dimension is the index of the evaluation point. 1399 1400 The ordering of the basis functions is not graded, so the basis functions are not nested by degree like `PetscDTPKDEvalJet()`. 1401 The basis functions are not an L2-orthonormal basis on any particular domain. 1402 1403 The implementation is based on the description of the trimmed polynomials up to degree r as 1404 the direct sum of polynomials up to degree (r-1) and the Koszul differential applied to 1405 homogeneous polynomials of degree (r-1). 1406 1407 .seealso: `PetscDTPKDEvalJet()`, `PetscDTPTrimmedSize()` 1408 @*/ 1409 PetscErrorCode PetscDTPTrimmedEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[]) 1410 { 1411 PetscFunctionBegin; 1412 PetscCall(PetscDTPTrimmedEvalJet_Internal(dim, npoints, points, degree, formDegree, jetDegree, p)); 1413 PetscFunctionReturn(PETSC_SUCCESS); 1414 } 1415 1416 /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V 1417 * with lds n; diag and subdiag are overwritten */ 1418 static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[], PetscReal eigs[], PetscScalar V[]) 1419 { 1420 char jobz = 'V'; /* eigenvalues and eigenvectors */ 1421 char range = 'A'; /* all eigenvalues will be found */ 1422 PetscReal VL = 0.; /* ignored because range is 'A' */ 1423 PetscReal VU = 0.; /* ignored because range is 'A' */ 1424 PetscBLASInt IL = 0; /* ignored because range is 'A' */ 1425 PetscBLASInt IU = 0; /* ignored because range is 'A' */ 1426 PetscReal abstol = 0.; /* unused */ 1427 PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */ 1428 PetscBLASInt *isuppz; 1429 PetscBLASInt lwork, liwork; 1430 PetscReal workquery; 1431 PetscBLASInt iworkquery; 1432 PetscBLASInt *iwork; 1433 PetscBLASInt info; 1434 PetscReal *work = NULL; 1435 1436 PetscFunctionBegin; 1437 #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG) 1438 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found"); 1439 #endif 1440 PetscCall(PetscBLASIntCast(n, &bn)); 1441 PetscCall(PetscBLASIntCast(n, &ldz)); 1442 #if !defined(PETSC_MISSING_LAPACK_STEGR) 1443 PetscCall(PetscMalloc1(2 * n, &isuppz)); 1444 lwork = -1; 1445 liwork = -1; 1446 PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, &workquery, &lwork, &iworkquery, &liwork, &info)); 1447 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEGR error"); 1448 lwork = (PetscBLASInt)workquery; 1449 liwork = iworkquery; 1450 PetscCall(PetscMalloc2(lwork, &work, liwork, &iwork)); 1451 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 1452 PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, work, &lwork, iwork, &liwork, &info)); 1453 PetscCall(PetscFPTrapPop()); 1454 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEGR error"); 1455 PetscCall(PetscFree2(work, iwork)); 1456 PetscCall(PetscFree(isuppz)); 1457 #elif !defined(PETSC_MISSING_LAPACK_STEQR) 1458 jobz = 'I'; /* Compute eigenvalues and eigenvectors of the 1459 tridiagonal matrix. Z is initialized to the identity 1460 matrix. */ 1461 PetscCall(PetscMalloc1(PetscMax(1, 2 * n - 2), &work)); 1462 PetscCallBLAS("LAPACKsteqr", LAPACKsteqr_("I", &bn, diag, subdiag, V, &ldz, work, &info)); 1463 PetscCall(PetscFPTrapPop()); 1464 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEQR error"); 1465 PetscCall(PetscFree(work)); 1466 PetscCall(PetscArraycpy(eigs, diag, n)); 1467 #endif 1468 PetscFunctionReturn(PETSC_SUCCESS); 1469 } 1470 1471 /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi 1472 * quadrature rules on the interval [-1, 1] */ 1473 static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw) 1474 { 1475 PetscReal twoab1; 1476 PetscInt m = n - 2; 1477 PetscReal a = alpha + 1.; 1478 PetscReal b = beta + 1.; 1479 PetscReal gra, grb; 1480 1481 PetscFunctionBegin; 1482 twoab1 = PetscPowReal(2., a + b - 1.); 1483 #if defined(PETSC_HAVE_LGAMMA) 1484 grb = PetscExpReal(2. * PetscLGamma(b + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + a + 1.) - (PetscLGamma(m + b + 1) + PetscLGamma(m + a + b + 1.))); 1485 gra = PetscExpReal(2. * PetscLGamma(a + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + b + 1.) - (PetscLGamma(m + a + 1) + PetscLGamma(m + a + b + 1.))); 1486 #else 1487 { 1488 PetscInt alphai = (PetscInt)alpha; 1489 PetscInt betai = (PetscInt)beta; 1490 1491 if ((PetscReal)alphai == alpha && (PetscReal)betai == beta) { 1492 PetscReal binom1, binom2; 1493 1494 PetscCall(PetscDTBinomial(m + b, b, &binom1)); 1495 PetscCall(PetscDTBinomial(m + a + b, b, &binom2)); 1496 grb = 1. / (binom1 * binom2); 1497 PetscCall(PetscDTBinomial(m + a, a, &binom1)); 1498 PetscCall(PetscDTBinomial(m + a + b, a, &binom2)); 1499 gra = 1. / (binom1 * binom2); 1500 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable."); 1501 } 1502 #endif 1503 *leftw = twoab1 * grb / b; 1504 *rightw = twoab1 * gra / a; 1505 PetscFunctionReturn(PETSC_SUCCESS); 1506 } 1507 1508 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x. 1509 Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */ 1510 static inline PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 1511 { 1512 PetscReal pn1, pn2; 1513 PetscReal cnm1, cnm1x, cnm2; 1514 PetscInt k; 1515 1516 PetscFunctionBegin; 1517 if (!n) { 1518 *P = 1.0; 1519 PetscFunctionReturn(PETSC_SUCCESS); 1520 } 1521 PetscDTJacobiRecurrence_Internal(1, a, b, cnm1, cnm1x, cnm2); 1522 pn2 = 1.; 1523 pn1 = cnm1 + cnm1x * x; 1524 if (n == 1) { 1525 *P = pn1; 1526 PetscFunctionReturn(PETSC_SUCCESS); 1527 } 1528 *P = 0.0; 1529 for (k = 2; k < n + 1; ++k) { 1530 PetscDTJacobiRecurrence_Internal(k, a, b, cnm1, cnm1x, cnm2); 1531 1532 *P = (cnm1 + cnm1x * x) * pn1 - cnm2 * pn2; 1533 pn2 = pn1; 1534 pn1 = *P; 1535 } 1536 PetscFunctionReturn(PETSC_SUCCESS); 1537 } 1538 1539 /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */ 1540 static inline PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P) 1541 { 1542 PetscReal nP; 1543 PetscInt i; 1544 1545 PetscFunctionBegin; 1546 *P = 0.0; 1547 if (k > n) PetscFunctionReturn(PETSC_SUCCESS); 1548 PetscCall(PetscDTComputeJacobi(a + k, b + k, n - k, x, &nP)); 1549 for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5; 1550 *P = nP; 1551 PetscFunctionReturn(PETSC_SUCCESS); 1552 } 1553 1554 static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[]) 1555 { 1556 PetscInt maxIter = 100; 1557 PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON)); 1558 PetscReal a1, a6, gf; 1559 PetscInt k; 1560 1561 PetscFunctionBegin; 1562 a1 = PetscPowReal(2.0, a + b + 1); 1563 #if defined(PETSC_HAVE_LGAMMA) 1564 { 1565 PetscReal a2, a3, a4, a5; 1566 a2 = PetscLGamma(a + npoints + 1); 1567 a3 = PetscLGamma(b + npoints + 1); 1568 a4 = PetscLGamma(a + b + npoints + 1); 1569 a5 = PetscLGamma(npoints + 1); 1570 gf = PetscExpReal(a2 + a3 - (a4 + a5)); 1571 } 1572 #else 1573 { 1574 PetscInt ia, ib; 1575 1576 ia = (PetscInt)a; 1577 ib = (PetscInt)b; 1578 gf = 1.; 1579 if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */ 1580 for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k); 1581 } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */ 1582 for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k); 1583 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable."); 1584 } 1585 #endif 1586 1587 a6 = a1 * gf; 1588 /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses. 1589 Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */ 1590 for (k = 0; k < npoints; ++k) { 1591 PetscReal r = PetscCosReal(PETSC_PI * (1. - (4. * k + 3. + 2. * b) / (4. * npoints + 2. * (a + b + 1.)))), dP; 1592 PetscInt j; 1593 1594 if (k > 0) r = 0.5 * (r + x[k - 1]); 1595 for (j = 0; j < maxIter; ++j) { 1596 PetscReal s = 0.0, delta, f, fp; 1597 PetscInt i; 1598 1599 for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]); 1600 PetscCall(PetscDTComputeJacobi(a, b, npoints, r, &f)); 1601 PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp)); 1602 delta = f / (fp - f * s); 1603 r = r - delta; 1604 if (PetscAbsReal(delta) < eps) break; 1605 } 1606 x[k] = r; 1607 PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP)); 1608 w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP); 1609 } 1610 PetscFunctionReturn(PETSC_SUCCESS); 1611 } 1612 1613 /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi 1614 * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */ 1615 static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s) 1616 { 1617 PetscInt i; 1618 1619 PetscFunctionBegin; 1620 for (i = 0; i < nPoints; i++) { 1621 PetscReal A, B, C; 1622 1623 PetscDTJacobiRecurrence_Internal(i + 1, a, b, A, B, C); 1624 d[i] = -A / B; 1625 if (i) s[i - 1] *= C / B; 1626 if (i < nPoints - 1) s[i] = 1. / B; 1627 } 1628 PetscFunctionReturn(PETSC_SUCCESS); 1629 } 1630 1631 static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w) 1632 { 1633 PetscReal mu0; 1634 PetscReal ga, gb, gab; 1635 PetscInt i; 1636 1637 PetscFunctionBegin; 1638 PetscCall(PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite)); 1639 1640 #if defined(PETSC_HAVE_TGAMMA) 1641 ga = PetscTGamma(a + 1); 1642 gb = PetscTGamma(b + 1); 1643 gab = PetscTGamma(a + b + 2); 1644 #else 1645 { 1646 PetscInt ia, ib; 1647 1648 ia = (PetscInt)a; 1649 ib = (PetscInt)b; 1650 if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */ 1651 PetscCall(PetscDTFactorial(ia, &ga)); 1652 PetscCall(PetscDTFactorial(ib, &gb)); 1653 PetscCall(PetscDTFactorial(ia + ib + 1, &gab)); 1654 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "tgamma() - math routine is unavailable."); 1655 } 1656 #endif 1657 mu0 = PetscPowReal(2., a + b + 1.) * ga * gb / gab; 1658 1659 #if defined(PETSCDTGAUSSIANQUADRATURE_EIG) 1660 { 1661 PetscReal *diag, *subdiag; 1662 PetscScalar *V; 1663 1664 PetscCall(PetscMalloc2(npoints, &diag, npoints, &subdiag)); 1665 PetscCall(PetscMalloc1(npoints * npoints, &V)); 1666 PetscCall(PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag)); 1667 for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]); 1668 PetscCall(PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V)); 1669 for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0; 1670 PetscCall(PetscFree(V)); 1671 PetscCall(PetscFree2(diag, subdiag)); 1672 } 1673 #else 1674 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found"); 1675 #endif 1676 { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the 1677 eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that 1678 the eigenvalues are sorted */ 1679 PetscBool sorted; 1680 1681 PetscCall(PetscSortedReal(npoints, x, &sorted)); 1682 if (!sorted) { 1683 PetscInt *order, i; 1684 PetscReal *tmp; 1685 1686 PetscCall(PetscMalloc2(npoints, &order, npoints, &tmp)); 1687 for (i = 0; i < npoints; i++) order[i] = i; 1688 PetscCall(PetscSortRealWithPermutation(npoints, x, order)); 1689 PetscCall(PetscArraycpy(tmp, x, npoints)); 1690 for (i = 0; i < npoints; i++) x[i] = tmp[order[i]]; 1691 PetscCall(PetscArraycpy(tmp, w, npoints)); 1692 for (i = 0; i < npoints; i++) w[i] = tmp[order[i]]; 1693 PetscCall(PetscFree2(order, tmp)); 1694 } 1695 } 1696 PetscFunctionReturn(PETSC_SUCCESS); 1697 } 1698 1699 static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1700 { 1701 PetscFunctionBegin; 1702 PetscCheck(npoints >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive"); 1703 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1704 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1."); 1705 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1."); 1706 1707 if (newton) PetscCall(PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w)); 1708 else PetscCall(PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w)); 1709 if (alpha == beta) { /* symmetrize */ 1710 PetscInt i; 1711 for (i = 0; i < (npoints + 1) / 2; i++) { 1712 PetscInt j = npoints - 1 - i; 1713 PetscReal xi = x[i]; 1714 PetscReal xj = x[j]; 1715 PetscReal wi = w[i]; 1716 PetscReal wj = w[j]; 1717 1718 x[i] = (xi - xj) / 2.; 1719 x[j] = (xj - xi) / 2.; 1720 w[i] = w[j] = (wi + wj) / 2.; 1721 } 1722 } 1723 PetscFunctionReturn(PETSC_SUCCESS); 1724 } 1725 1726 /*@ 1727 PetscDTGaussJacobiQuadrature - quadrature for the interval $[a, b]$ with the weight function 1728 $(x-a)^\alpha (x-b)^\beta$. 1729 1730 Not Collective 1731 1732 Input Parameters: 1733 + npoints - the number of points in the quadrature rule 1734 . a - the left endpoint of the interval 1735 . b - the right endpoint of the interval 1736 . alpha - the left exponent 1737 - beta - the right exponent 1738 1739 Output Parameters: 1740 + x - array of length `npoints`, the locations of the quadrature points 1741 - w - array of length `npoints`, the weights of the quadrature points 1742 1743 Level: intermediate 1744 1745 Note: 1746 This quadrature rule is exact for polynomials up to degree 2*`npoints` - 1. 1747 1748 .seealso: `PetscDTGaussQuadrature()` 1749 @*/ 1750 PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1751 { 1752 PetscInt i; 1753 1754 PetscFunctionBegin; 1755 PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal)); 1756 if (a != -1. || b != 1.) { /* shift */ 1757 for (i = 0; i < npoints; i++) { 1758 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1759 w[i] *= (b - a) / 2.; 1760 } 1761 } 1762 PetscFunctionReturn(PETSC_SUCCESS); 1763 } 1764 1765 static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1766 { 1767 PetscInt i; 1768 1769 PetscFunctionBegin; 1770 PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive"); 1771 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1772 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1."); 1773 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1."); 1774 1775 x[0] = -1.; 1776 x[npoints - 1] = 1.; 1777 if (npoints > 2) PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints - 2, alpha + 1., beta + 1., &x[1], &w[1], newton)); 1778 for (i = 1; i < npoints - 1; i++) w[i] /= (1. - x[i] * x[i]); 1779 PetscCall(PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints - 1])); 1780 PetscFunctionReturn(PETSC_SUCCESS); 1781 } 1782 1783 /*@ 1784 PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval $[a, b]$ with the weight function 1785 $(x-a)^\alpha (x-b)^\beta$, with endpoints `a` and `b` included as quadrature points. 1786 1787 Not Collective 1788 1789 Input Parameters: 1790 + npoints - the number of points in the quadrature rule 1791 . a - the left endpoint of the interval 1792 . b - the right endpoint of the interval 1793 . alpha - the left exponent 1794 - beta - the right exponent 1795 1796 Output Parameters: 1797 + x - array of length `npoints`, the locations of the quadrature points 1798 - w - array of length `npoints`, the weights of the quadrature points 1799 1800 Level: intermediate 1801 1802 Note: 1803 This quadrature rule is exact for polynomials up to degree 2*`npoints` - 3. 1804 1805 .seealso: `PetscDTGaussJacobiQuadrature()` 1806 @*/ 1807 PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1808 { 1809 PetscInt i; 1810 1811 PetscFunctionBegin; 1812 PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal)); 1813 if (a != -1. || b != 1.) { /* shift */ 1814 for (i = 0; i < npoints; i++) { 1815 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1816 w[i] *= (b - a) / 2.; 1817 } 1818 } 1819 PetscFunctionReturn(PETSC_SUCCESS); 1820 } 1821 1822 /*@ 1823 PetscDTGaussQuadrature - create Gauss-Legendre quadrature 1824 1825 Not Collective 1826 1827 Input Parameters: 1828 + npoints - number of points 1829 . a - left end of interval (often-1) 1830 - b - right end of interval (often +1) 1831 1832 Output Parameters: 1833 + x - quadrature points 1834 - w - quadrature weights 1835 1836 Level: intermediate 1837 1838 Note: 1839 See {cite}`golub1969calculation` 1840 1841 .seealso: `PetscDTLegendreEval()`, `PetscDTGaussJacobiQuadrature()` 1842 @*/ 1843 PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[]) 1844 { 1845 PetscInt i; 1846 1847 PetscFunctionBegin; 1848 PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal)); 1849 if (a != -1. || b != 1.) { /* shift */ 1850 for (i = 0; i < npoints; i++) { 1851 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1852 w[i] *= (b - a) / 2.; 1853 } 1854 } 1855 PetscFunctionReturn(PETSC_SUCCESS); 1856 } 1857 1858 /*@ 1859 PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre 1860 nodes of a given size on the domain $[-1,1]$ 1861 1862 Not Collective 1863 1864 Input Parameters: 1865 + npoints - number of grid nodes 1866 - type - `PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA` or `PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON` 1867 1868 Output Parameters: 1869 + x - quadrature points, pass in an array of length `npoints` 1870 - w - quadrature weights, pass in an array of length `npoints` 1871 1872 Level: intermediate 1873 1874 Notes: 1875 For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not 1876 close enough to the desired solution 1877 1878 These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes 1879 1880 See <https://epubs.siam.org/doi/abs/10.1137/110855442> <https://epubs.siam.org/doi/abs/10.1137/120889873> for better ways to compute GLL nodes 1881 1882 .seealso: `PetscDTGaussQuadrature()`, `PetscGaussLobattoLegendreCreateType` 1883 1884 @*/ 1885 PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints, PetscGaussLobattoLegendreCreateType type, PetscReal x[], PetscReal w[]) 1886 { 1887 PetscBool newton; 1888 1889 PetscFunctionBegin; 1890 PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide at least 2 grid points per element"); 1891 newton = (PetscBool)(type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON); 1892 PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton)); 1893 PetscFunctionReturn(PETSC_SUCCESS); 1894 } 1895 1896 /*@ 1897 PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature 1898 1899 Not Collective 1900 1901 Input Parameters: 1902 + dim - The spatial dimension 1903 . Nc - The number of components 1904 . npoints - number of points in one dimension 1905 . a - left end of interval (often-1) 1906 - b - right end of interval (often +1) 1907 1908 Output Parameter: 1909 . q - A `PetscQuadrature` object 1910 1911 Level: intermediate 1912 1913 .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()` 1914 @*/ 1915 PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 1916 { 1917 DMPolytopeType ct; 1918 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints * PetscSqr(npoints) : PetscSqr(npoints) : npoints; 1919 PetscReal *x, *w, *xw, *ww; 1920 1921 PetscFunctionBegin; 1922 PetscCall(PetscMalloc1(totpoints * dim, &x)); 1923 PetscCall(PetscMalloc1(totpoints * Nc, &w)); 1924 /* Set up the Golub-Welsch system */ 1925 switch (dim) { 1926 case 0: 1927 ct = DM_POLYTOPE_POINT; 1928 PetscCall(PetscFree(x)); 1929 PetscCall(PetscFree(w)); 1930 PetscCall(PetscMalloc1(1, &x)); 1931 PetscCall(PetscMalloc1(Nc, &w)); 1932 totpoints = 1; 1933 x[0] = 0.0; 1934 for (PetscInt c = 0; c < Nc; ++c) w[c] = 1.0; 1935 break; 1936 case 1: 1937 ct = DM_POLYTOPE_SEGMENT; 1938 PetscCall(PetscMalloc1(npoints, &ww)); 1939 PetscCall(PetscDTGaussQuadrature(npoints, a, b, x, ww)); 1940 for (PetscInt i = 0; i < npoints; ++i) 1941 for (PetscInt c = 0; c < Nc; ++c) w[i * Nc + c] = ww[i]; 1942 PetscCall(PetscFree(ww)); 1943 break; 1944 case 2: 1945 ct = DM_POLYTOPE_QUADRILATERAL; 1946 PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww)); 1947 PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww)); 1948 for (PetscInt i = 0; i < npoints; ++i) { 1949 for (PetscInt j = 0; j < npoints; ++j) { 1950 x[(i * npoints + j) * dim + 0] = xw[i]; 1951 x[(i * npoints + j) * dim + 1] = xw[j]; 1952 for (PetscInt c = 0; c < Nc; ++c) w[(i * npoints + j) * Nc + c] = ww[i] * ww[j]; 1953 } 1954 } 1955 PetscCall(PetscFree2(xw, ww)); 1956 break; 1957 case 3: 1958 ct = DM_POLYTOPE_HEXAHEDRON; 1959 PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww)); 1960 PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww)); 1961 for (PetscInt i = 0; i < npoints; ++i) { 1962 for (PetscInt j = 0; j < npoints; ++j) { 1963 for (PetscInt k = 0; k < npoints; ++k) { 1964 x[((i * npoints + j) * npoints + k) * dim + 0] = xw[i]; 1965 x[((i * npoints + j) * npoints + k) * dim + 1] = xw[j]; 1966 x[((i * npoints + j) * npoints + k) * dim + 2] = xw[k]; 1967 for (PetscInt c = 0; c < Nc; ++c) w[((i * npoints + j) * npoints + k) * Nc + c] = ww[i] * ww[j] * ww[k]; 1968 } 1969 } 1970 } 1971 PetscCall(PetscFree2(xw, ww)); 1972 break; 1973 default: 1974 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %" PetscInt_FMT, dim); 1975 } 1976 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 1977 PetscCall(PetscQuadratureSetCellType(*q, ct)); 1978 PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1)); 1979 PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w)); 1980 PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "GaussTensor")); 1981 PetscFunctionReturn(PETSC_SUCCESS); 1982 } 1983 1984 /*@ 1985 PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex {cite}`karniadakis2005spectral` 1986 1987 Not Collective 1988 1989 Input Parameters: 1990 + dim - The simplex dimension 1991 . Nc - The number of components 1992 . npoints - The number of points in one dimension 1993 . a - left end of interval (often-1) 1994 - b - right end of interval (often +1) 1995 1996 Output Parameter: 1997 . q - A `PetscQuadrature` object 1998 1999 Level: intermediate 2000 2001 Note: 2002 For `dim` == 1, this is Gauss-Legendre quadrature 2003 2004 .seealso: `PetscDTGaussTensorQuadrature()`, `PetscDTGaussQuadrature()` 2005 @*/ 2006 PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 2007 { 2008 DMPolytopeType ct; 2009 PetscInt totpoints; 2010 PetscReal *p1, *w1; 2011 PetscReal *x, *w; 2012 2013 PetscFunctionBegin; 2014 PetscCheck(!(a != -1.0) && !(b != 1.0), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now"); 2015 switch (dim) { 2016 case 0: 2017 ct = DM_POLYTOPE_POINT; 2018 break; 2019 case 1: 2020 ct = DM_POLYTOPE_SEGMENT; 2021 break; 2022 case 2: 2023 ct = DM_POLYTOPE_TRIANGLE; 2024 break; 2025 case 3: 2026 ct = DM_POLYTOPE_TETRAHEDRON; 2027 break; 2028 default: 2029 ct = DM_POLYTOPE_UNKNOWN; 2030 } 2031 totpoints = 1; 2032 for (PetscInt i = 0; i < dim; ++i) totpoints *= npoints; 2033 PetscCall(PetscMalloc1(totpoints * dim, &x)); 2034 PetscCall(PetscMalloc1(totpoints * Nc, &w)); 2035 PetscCall(PetscMalloc2(npoints, &p1, npoints, &w1)); 2036 for (PetscInt i = 0; i < totpoints * Nc; ++i) w[i] = 1.; 2037 for (PetscInt i = 0, totprev = 1, totrem = totpoints / npoints; i < dim; ++i) { 2038 PetscReal mul; 2039 2040 mul = PetscPowReal(2., -i); 2041 PetscCall(PetscDTGaussJacobiQuadrature(npoints, -1., 1., i, 0.0, p1, w1)); 2042 for (PetscInt pt = 0, l = 0; l < totprev; l++) { 2043 for (PetscInt j = 0; j < npoints; j++) { 2044 for (PetscInt m = 0; m < totrem; m++, pt++) { 2045 for (PetscInt k = 0; k < i; k++) x[pt * dim + k] = (x[pt * dim + k] + 1.) * (1. - p1[j]) * 0.5 - 1.; 2046 x[pt * dim + i] = p1[j]; 2047 for (PetscInt c = 0; c < Nc; c++) w[pt * Nc + c] *= mul * w1[j]; 2048 } 2049 } 2050 } 2051 totprev *= npoints; 2052 totrem /= npoints; 2053 } 2054 PetscCall(PetscFree2(p1, w1)); 2055 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 2056 PetscCall(PetscQuadratureSetCellType(*q, ct)); 2057 PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1)); 2058 PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w)); 2059 PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "StroudConical")); 2060 PetscFunctionReturn(PETSC_SUCCESS); 2061 } 2062 2063 static PetscBool MinSymTriQuadCite = PETSC_FALSE; 2064 const char MinSymTriQuadCitation[] = "@article{WitherdenVincent2015,\n" 2065 " title = {On the identification of symmetric quadrature rules for finite element methods},\n" 2066 " journal = {Computers & Mathematics with Applications},\n" 2067 " volume = {69},\n" 2068 " number = {10},\n" 2069 " pages = {1232-1241},\n" 2070 " year = {2015},\n" 2071 " issn = {0898-1221},\n" 2072 " doi = {10.1016/j.camwa.2015.03.017},\n" 2073 " url = {https://www.sciencedirect.com/science/article/pii/S0898122115001224},\n" 2074 " author = {F.D. Witherden and P.E. Vincent},\n" 2075 "}\n"; 2076 2077 #include "petscdttriquadrules.h" 2078 2079 static PetscBool MinSymTetQuadCite = PETSC_FALSE; 2080 const char MinSymTetQuadCitation[] = "@article{JaskowiecSukumar2021\n" 2081 " author = {Jaskowiec, Jan and Sukumar, N.},\n" 2082 " title = {High-order symmetric cubature rules for tetrahedra and pyramids},\n" 2083 " journal = {International Journal for Numerical Methods in Engineering},\n" 2084 " volume = {122},\n" 2085 " number = {1},\n" 2086 " pages = {148-171},\n" 2087 " doi = {10.1002/nme.6528},\n" 2088 " url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6528},\n" 2089 " eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6528},\n" 2090 " year = {2021}\n" 2091 "}\n"; 2092 2093 #include "petscdttetquadrules.h" 2094 2095 // https://en.wikipedia.org/wiki/Partition_(number_theory) 2096 static PetscErrorCode PetscDTPartitionNumber(PetscInt n, PetscInt *p) 2097 { 2098 // sequence A000041 in the OEIS 2099 const PetscInt partition[] = {1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135, 176, 231, 297, 385, 490, 627, 792, 1002, 1255, 1575, 1958, 2436, 3010, 3718, 4565, 5604}; 2100 PetscInt tabulated_max = PETSC_STATIC_ARRAY_LENGTH(partition) - 1; 2101 2102 PetscFunctionBegin; 2103 PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Partition number not defined for negative number %" PetscInt_FMT, n); 2104 // not implementing the pentagonal number recurrence, we don't need partition numbers for n that high 2105 PetscCheck(n <= tabulated_max, PETSC_COMM_SELF, PETSC_ERR_SUP, "Partition numbers only tabulated up to %" PetscInt_FMT ", not computed for %" PetscInt_FMT, tabulated_max, n); 2106 *p = partition[n]; 2107 PetscFunctionReturn(PETSC_SUCCESS); 2108 } 2109 2110 /*@ 2111 PetscDTSimplexQuadrature - Create a quadrature rule for a simplex that exactly integrates polynomials up to a given degree. 2112 2113 Not Collective 2114 2115 Input Parameters: 2116 + dim - The spatial dimension of the simplex (1 = segment, 2 = triangle, 3 = tetrahedron) 2117 . degree - The largest polynomial degree that is required to be integrated exactly 2118 - type - left end of interval (often-1) 2119 2120 Output Parameter: 2121 . quad - A `PetscQuadrature` object for integration over the biunit simplex 2122 (defined by the bounds $x_i >= -1$ and $\sum_i x_i <= 2 - d$) that is exact for 2123 polynomials up to the given degree 2124 2125 Level: intermediate 2126 2127 .seealso: `PetscDTSimplexQuadratureType`, `PetscDTGaussQuadrature()`, `PetscDTStroudCononicalQuadrature()`, `PetscQuadrature` 2128 @*/ 2129 PetscErrorCode PetscDTSimplexQuadrature(PetscInt dim, PetscInt degree, PetscDTSimplexQuadratureType type, PetscQuadrature *quad) 2130 { 2131 PetscDTSimplexQuadratureType orig_type = type; 2132 2133 PetscFunctionBegin; 2134 PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative dimension %" PetscInt_FMT, dim); 2135 PetscCheck(degree >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %" PetscInt_FMT, degree); 2136 if (type == PETSCDTSIMPLEXQUAD_DEFAULT) type = PETSCDTSIMPLEXQUAD_MINSYM; 2137 if (type == PETSCDTSIMPLEXQUAD_CONIC || dim < 2) { 2138 PetscInt points_per_dim = (degree + 2) / 2; // ceil((degree + 1) / 2); 2139 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, points_per_dim, -1, 1, quad)); 2140 } else { 2141 DMPolytopeType ct; 2142 PetscInt n = dim + 1; 2143 PetscInt fact = 1; 2144 PetscInt *part, *perm; 2145 PetscInt p = 0; 2146 PetscInt max_degree; 2147 const PetscInt *nodes_per_type = NULL; 2148 const PetscInt *all_num_full_nodes = NULL; 2149 const PetscReal **weights_list = NULL; 2150 const PetscReal **compact_nodes_list = NULL; 2151 const char *citation = NULL; 2152 PetscBool *cited = NULL; 2153 2154 switch (dim) { 2155 case 0: 2156 ct = DM_POLYTOPE_POINT; 2157 break; 2158 case 1: 2159 ct = DM_POLYTOPE_SEGMENT; 2160 break; 2161 case 2: 2162 ct = DM_POLYTOPE_TRIANGLE; 2163 break; 2164 case 3: 2165 ct = DM_POLYTOPE_TETRAHEDRON; 2166 break; 2167 default: 2168 ct = DM_POLYTOPE_UNKNOWN; 2169 } 2170 switch (dim) { 2171 case 2: 2172 cited = &MinSymTriQuadCite; 2173 citation = MinSymTriQuadCitation; 2174 max_degree = PetscDTWVTriQuad_max_degree; 2175 nodes_per_type = PetscDTWVTriQuad_num_orbits; 2176 all_num_full_nodes = PetscDTWVTriQuad_num_nodes; 2177 weights_list = PetscDTWVTriQuad_weights; 2178 compact_nodes_list = PetscDTWVTriQuad_orbits; 2179 break; 2180 case 3: 2181 cited = &MinSymTetQuadCite; 2182 citation = MinSymTetQuadCitation; 2183 max_degree = PetscDTJSTetQuad_max_degree; 2184 nodes_per_type = PetscDTJSTetQuad_num_orbits; 2185 all_num_full_nodes = PetscDTJSTetQuad_num_nodes; 2186 weights_list = PetscDTJSTetQuad_weights; 2187 compact_nodes_list = PetscDTJSTetQuad_orbits; 2188 break; 2189 default: 2190 max_degree = -1; 2191 break; 2192 } 2193 2194 if (degree > max_degree) { 2195 if (orig_type == PETSCDTSIMPLEXQUAD_DEFAULT) { 2196 // fall back to conic 2197 PetscCall(PetscDTSimplexQuadrature(dim, degree, PETSCDTSIMPLEXQUAD_CONIC, quad)); 2198 PetscFunctionReturn(PETSC_SUCCESS); 2199 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Minimal symmetric quadrature for dim %" PetscInt_FMT ", degree %" PetscInt_FMT " unsupported", dim, degree); 2200 } 2201 2202 PetscCall(PetscCitationsRegister(citation, cited)); 2203 2204 PetscCall(PetscDTPartitionNumber(n, &p)); 2205 for (PetscInt d = 2; d <= n; d++) fact *= d; 2206 2207 PetscInt num_full_nodes = all_num_full_nodes[degree]; 2208 const PetscReal *all_compact_nodes = compact_nodes_list[degree]; 2209 const PetscReal *all_compact_weights = weights_list[degree]; 2210 nodes_per_type = &nodes_per_type[p * degree]; 2211 2212 PetscReal *points; 2213 PetscReal *counts; 2214 PetscReal *weights; 2215 PetscReal *bary_to_biunit; // row-major transformation of barycentric coordinate to biunit 2216 PetscQuadrature q; 2217 2218 // compute the transformation 2219 PetscCall(PetscMalloc1(n * dim, &bary_to_biunit)); 2220 for (PetscInt d = 0; d < dim; d++) { 2221 for (PetscInt b = 0; b < n; b++) bary_to_biunit[d * n + b] = (d == b) ? 1.0 : -1.0; 2222 } 2223 2224 PetscCall(PetscMalloc3(n, &part, n, &perm, n, &counts)); 2225 PetscCall(PetscCalloc1(num_full_nodes * dim, &points)); 2226 PetscCall(PetscMalloc1(num_full_nodes, &weights)); 2227 2228 // (0, 0, ...) is the first partition lexicographically 2229 PetscCall(PetscArrayzero(part, n)); 2230 PetscCall(PetscArrayzero(counts, n)); 2231 counts[0] = n; 2232 2233 // for each partition 2234 for (PetscInt s = 0, node_offset = 0; s < p; s++) { 2235 PetscInt num_compact_coords = part[n - 1] + 1; 2236 2237 const PetscReal *compact_nodes = all_compact_nodes; 2238 const PetscReal *compact_weights = all_compact_weights; 2239 all_compact_nodes += num_compact_coords * nodes_per_type[s]; 2240 all_compact_weights += nodes_per_type[s]; 2241 2242 // for every permutation of the vertices 2243 for (PetscInt f = 0; f < fact; f++) { 2244 PetscCall(PetscDTEnumPerm(n, f, perm, NULL)); 2245 2246 // check if it is a valid permutation 2247 PetscInt digit; 2248 for (digit = 1; digit < n; digit++) { 2249 // skip permutations that would duplicate a node because it has a smaller symmetry group 2250 if (part[digit - 1] == part[digit] && perm[digit - 1] > perm[digit]) break; 2251 } 2252 if (digit < n) continue; 2253 2254 // create full nodes from this permutation of the compact nodes 2255 PetscReal *full_nodes = &points[node_offset * dim]; 2256 PetscReal *full_weights = &weights[node_offset]; 2257 2258 PetscCall(PetscArraycpy(full_weights, compact_weights, nodes_per_type[s])); 2259 for (PetscInt b = 0; b < n; b++) { 2260 for (PetscInt d = 0; d < dim; d++) { 2261 for (PetscInt node = 0; node < nodes_per_type[s]; node++) full_nodes[node * dim + d] += bary_to_biunit[d * n + perm[b]] * compact_nodes[node * num_compact_coords + part[b]]; 2262 } 2263 } 2264 node_offset += nodes_per_type[s]; 2265 } 2266 2267 if (s < p - 1) { // Generate the next partition 2268 /* A partition is described by the number of coordinates that are in 2269 * each set of duplicates (counts) and redundantly by mapping each 2270 * index to its set of duplicates (part) 2271 * 2272 * Counts should always be in nonincreasing order 2273 * 2274 * We want to generate the partitions lexically by part, which means 2275 * finding the last index where count > 1 and reducing by 1. 2276 * 2277 * For the new counts beyond that index, we eagerly assign the remaining 2278 * capacity of the partition to smaller indices (ensures lexical ordering), 2279 * while respecting the nonincreasing invariant of the counts 2280 */ 2281 PetscInt last_digit = part[n - 1]; 2282 PetscInt last_digit_with_extra = last_digit; 2283 while (counts[last_digit_with_extra] == 1) last_digit_with_extra--; 2284 PetscInt limit = --counts[last_digit_with_extra]; 2285 PetscInt total_to_distribute = last_digit - last_digit_with_extra + 1; 2286 for (PetscInt digit = last_digit_with_extra + 1; digit < n; digit++) { 2287 counts[digit] = PetscMin(limit, total_to_distribute); 2288 total_to_distribute -= PetscMin(limit, total_to_distribute); 2289 } 2290 for (PetscInt digit = 0, offset = 0; digit < n; digit++) { 2291 PetscInt count = counts[digit]; 2292 for (PetscInt c = 0; c < count; c++) part[offset++] = digit; 2293 } 2294 } 2295 } 2296 PetscCall(PetscFree3(part, perm, counts)); 2297 PetscCall(PetscFree(bary_to_biunit)); 2298 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &q)); 2299 PetscCall(PetscQuadratureSetCellType(q, ct)); 2300 PetscCall(PetscQuadratureSetOrder(q, degree)); 2301 PetscCall(PetscQuadratureSetData(q, dim, 1, num_full_nodes, points, weights)); 2302 *quad = q; 2303 } 2304 PetscFunctionReturn(PETSC_SUCCESS); 2305 } 2306 2307 /*@ 2308 PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell 2309 2310 Not Collective 2311 2312 Input Parameters: 2313 + dim - The cell dimension 2314 . level - The number of points in one dimension, $2^l$ 2315 . a - left end of interval (often-1) 2316 - b - right end of interval (often +1) 2317 2318 Output Parameter: 2319 . q - A `PetscQuadrature` object 2320 2321 Level: intermediate 2322 2323 .seealso: `PetscDTGaussTensorQuadrature()`, `PetscQuadrature` 2324 @*/ 2325 PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q) 2326 { 2327 DMPolytopeType ct; 2328 const PetscInt p = 16; /* Digits of precision in the evaluation */ 2329 const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */ 2330 const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */ 2331 const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */ 2332 PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */ 2333 PetscReal wk = 0.5 * PETSC_PI; /* Quadrature weight at x_k */ 2334 PetscReal *x, *w; 2335 PetscInt K, k, npoints; 2336 2337 PetscFunctionBegin; 2338 PetscCheck(dim <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not yet implemented", dim); 2339 PetscCheck(level, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits"); 2340 switch (dim) { 2341 case 0: 2342 ct = DM_POLYTOPE_POINT; 2343 break; 2344 case 1: 2345 ct = DM_POLYTOPE_SEGMENT; 2346 break; 2347 case 2: 2348 ct = DM_POLYTOPE_QUADRILATERAL; 2349 break; 2350 case 3: 2351 ct = DM_POLYTOPE_HEXAHEDRON; 2352 break; 2353 default: 2354 ct = DM_POLYTOPE_UNKNOWN; 2355 } 2356 /* Find K such that the weights are < 32 digits of precision */ 2357 for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2 * p; ++K) wk = 0.5 * h * PETSC_PI * PetscCoshReal(K * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(K * h))); 2358 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 2359 PetscCall(PetscQuadratureSetCellType(*q, ct)); 2360 PetscCall(PetscQuadratureSetOrder(*q, 2 * K + 1)); 2361 npoints = 2 * K - 1; 2362 PetscCall(PetscMalloc1(npoints * dim, &x)); 2363 PetscCall(PetscMalloc1(npoints, &w)); 2364 /* Center term */ 2365 x[0] = beta; 2366 w[0] = 0.5 * alpha * PETSC_PI; 2367 for (k = 1; k < K; ++k) { 2368 wk = 0.5 * alpha * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2369 xk = PetscTanhReal(0.5 * PETSC_PI * PetscSinhReal(k * h)); 2370 x[2 * k - 1] = -alpha * xk + beta; 2371 w[2 * k - 1] = wk; 2372 x[2 * k + 0] = alpha * xk + beta; 2373 w[2 * k + 0] = wk; 2374 } 2375 PetscCall(PetscQuadratureSetData(*q, dim, 1, npoints, x, w)); 2376 PetscFunctionReturn(PETSC_SUCCESS); 2377 } 2378 2379 PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2380 { 2381 const PetscInt p = 16; /* Digits of precision in the evaluation */ 2382 const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */ 2383 const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */ 2384 PetscReal h = 1.0; /* Step size, length between x_k */ 2385 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 2386 PetscReal osum = 0.0; /* Integral on last level */ 2387 PetscReal psum = 0.0; /* Integral on the level before the last level */ 2388 PetscReal sum; /* Integral on current level */ 2389 PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 2390 PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 2391 PetscReal wk; /* Quadrature weight at x_k */ 2392 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 2393 PetscInt d; /* Digits of precision in the integral */ 2394 2395 PetscFunctionBegin; 2396 PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 2397 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 2398 /* Center term */ 2399 func(&beta, ctx, &lval); 2400 sum = 0.5 * alpha * PETSC_PI * lval; 2401 /* */ 2402 do { 2403 PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4; 2404 PetscInt k = 1; 2405 2406 ++l; 2407 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */ 2408 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 2409 psum = osum; 2410 osum = sum; 2411 h *= 0.5; 2412 sum *= 0.5; 2413 do { 2414 wk = 0.5 * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2415 yk = 1.0 / (PetscExpReal(0.5 * PETSC_PI * PetscSinhReal(k * h)) * PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2416 lx = -alpha * (1.0 - yk) + beta; 2417 rx = alpha * (1.0 - yk) + beta; 2418 func(&lx, ctx, &lval); 2419 func(&rx, ctx, &rval); 2420 lterm = alpha * wk * lval; 2421 maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm); 2422 sum += lterm; 2423 rterm = alpha * wk * rval; 2424 maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm); 2425 sum += rterm; 2426 ++k; 2427 /* Only need to evaluate every other point on refined levels */ 2428 if (l != 1) ++k; 2429 } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */ 2430 2431 d1 = PetscLog10Real(PetscAbsReal(sum - osum)); 2432 d2 = PetscLog10Real(PetscAbsReal(sum - psum)); 2433 d3 = PetscLog10Real(maxTerm) - p; 2434 if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0; 2435 else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm))); 2436 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4))); 2437 } while (d < digits && l < 12); 2438 *sol = sum; 2439 PetscCall(PetscFPTrapPop()); 2440 PetscFunctionReturn(PETSC_SUCCESS); 2441 } 2442 2443 #if defined(PETSC_HAVE_MPFR) 2444 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2445 { 2446 const PetscInt safetyFactor = 2; /* Calculate abscissa until 2*p digits */ 2447 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 2448 mpfr_t alpha; /* Half-width of the integration interval */ 2449 mpfr_t beta; /* Center of the integration interval */ 2450 mpfr_t h; /* Step size, length between x_k */ 2451 mpfr_t osum; /* Integral on last level */ 2452 mpfr_t psum; /* Integral on the level before the last level */ 2453 mpfr_t sum; /* Integral on current level */ 2454 mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 2455 mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 2456 mpfr_t wk; /* Quadrature weight at x_k */ 2457 PetscReal lval, rval, rtmp; /* Terms in the quadature sum to the left and right of 0 */ 2458 PetscInt d; /* Digits of precision in the integral */ 2459 mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp; 2460 2461 PetscFunctionBegin; 2462 PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 2463 /* Create high precision storage */ 2464 mpfr_inits2(PetscCeilReal(safetyFactor * digits * PetscLogReal(10.) / PetscLogReal(2.)), alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL); 2465 /* Initialization */ 2466 mpfr_set_d(alpha, 0.5 * (b - a), MPFR_RNDN); 2467 mpfr_set_d(beta, 0.5 * (b + a), MPFR_RNDN); 2468 mpfr_set_d(osum, 0.0, MPFR_RNDN); 2469 mpfr_set_d(psum, 0.0, MPFR_RNDN); 2470 mpfr_set_d(h, 1.0, MPFR_RNDN); 2471 mpfr_const_pi(pi2, MPFR_RNDN); 2472 mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN); 2473 /* Center term */ 2474 rtmp = 0.5 * (b + a); 2475 func(&rtmp, ctx, &lval); 2476 mpfr_set(sum, pi2, MPFR_RNDN); 2477 mpfr_mul(sum, sum, alpha, MPFR_RNDN); 2478 mpfr_mul_d(sum, sum, lval, MPFR_RNDN); 2479 /* */ 2480 do { 2481 PetscReal d1, d2, d3, d4; 2482 PetscInt k = 1; 2483 2484 ++l; 2485 mpfr_set_d(maxTerm, 0.0, MPFR_RNDN); 2486 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */ 2487 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 2488 mpfr_set(psum, osum, MPFR_RNDN); 2489 mpfr_set(osum, sum, MPFR_RNDN); 2490 mpfr_mul_d(h, h, 0.5, MPFR_RNDN); 2491 mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN); 2492 do { 2493 mpfr_set_si(kh, k, MPFR_RNDN); 2494 mpfr_mul(kh, kh, h, MPFR_RNDN); 2495 /* Weight */ 2496 mpfr_set(wk, h, MPFR_RNDN); 2497 mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN); 2498 mpfr_mul(msinh, msinh, pi2, MPFR_RNDN); 2499 mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN); 2500 mpfr_cosh(tmp, msinh, MPFR_RNDN); 2501 mpfr_sqr(tmp, tmp, MPFR_RNDN); 2502 mpfr_mul(wk, wk, mcosh, MPFR_RNDN); 2503 mpfr_div(wk, wk, tmp, MPFR_RNDN); 2504 /* Abscissa */ 2505 mpfr_set_d(yk, 1.0, MPFR_RNDZ); 2506 mpfr_cosh(tmp, msinh, MPFR_RNDN); 2507 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 2508 mpfr_exp(tmp, msinh, MPFR_RNDN); 2509 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 2510 /* Quadrature points */ 2511 mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ); 2512 mpfr_mul(lx, lx, alpha, MPFR_RNDU); 2513 mpfr_add(lx, lx, beta, MPFR_RNDU); 2514 mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ); 2515 mpfr_mul(rx, rx, alpha, MPFR_RNDD); 2516 mpfr_add(rx, rx, beta, MPFR_RNDD); 2517 /* Evaluation */ 2518 rtmp = mpfr_get_d(lx, MPFR_RNDU); 2519 func(&rtmp, ctx, &lval); 2520 rtmp = mpfr_get_d(rx, MPFR_RNDD); 2521 func(&rtmp, ctx, &rval); 2522 /* Update */ 2523 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 2524 mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN); 2525 mpfr_add(sum, sum, tmp, MPFR_RNDN); 2526 mpfr_abs(tmp, tmp, MPFR_RNDN); 2527 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 2528 mpfr_set(curTerm, tmp, MPFR_RNDN); 2529 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 2530 mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN); 2531 mpfr_add(sum, sum, tmp, MPFR_RNDN); 2532 mpfr_abs(tmp, tmp, MPFR_RNDN); 2533 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 2534 mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN); 2535 ++k; 2536 /* Only need to evaluate every other point on refined levels */ 2537 if (l != 1) ++k; 2538 mpfr_log10(tmp, wk, MPFR_RNDN); 2539 mpfr_abs(tmp, tmp, MPFR_RNDN); 2540 } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor * digits); /* Only need to evaluate sum until weights are < 32 digits of precision */ 2541 mpfr_sub(tmp, sum, osum, MPFR_RNDN); 2542 mpfr_abs(tmp, tmp, MPFR_RNDN); 2543 mpfr_log10(tmp, tmp, MPFR_RNDN); 2544 d1 = mpfr_get_d(tmp, MPFR_RNDN); 2545 mpfr_sub(tmp, sum, psum, MPFR_RNDN); 2546 mpfr_abs(tmp, tmp, MPFR_RNDN); 2547 mpfr_log10(tmp, tmp, MPFR_RNDN); 2548 d2 = mpfr_get_d(tmp, MPFR_RNDN); 2549 mpfr_log10(tmp, maxTerm, MPFR_RNDN); 2550 d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits; 2551 mpfr_log10(tmp, curTerm, MPFR_RNDN); 2552 d4 = mpfr_get_d(tmp, MPFR_RNDN); 2553 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4))); 2554 } while (d < digits && l < 8); 2555 *sol = mpfr_get_d(sum, MPFR_RNDN); 2556 /* Cleanup */ 2557 mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL); 2558 PetscFunctionReturn(PETSC_SUCCESS); 2559 } 2560 #else 2561 2562 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2563 { 2564 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp"); 2565 } 2566 #endif 2567 2568 /*@ 2569 PetscDTTensorQuadratureCreate - create the tensor product quadrature from two lower-dimensional quadratures 2570 2571 Not Collective 2572 2573 Input Parameters: 2574 + q1 - The first quadrature 2575 - q2 - The second quadrature 2576 2577 Output Parameter: 2578 . q - A `PetscQuadrature` object 2579 2580 Level: intermediate 2581 2582 .seealso: `PetscQuadrature`, `PetscDTGaussTensorQuadrature()` 2583 @*/ 2584 PetscErrorCode PetscDTTensorQuadratureCreate(PetscQuadrature q1, PetscQuadrature q2, PetscQuadrature *q) 2585 { 2586 DMPolytopeType ct1, ct2, ct; 2587 const PetscReal *x1, *w1, *x2, *w2; 2588 PetscReal *x, *w; 2589 PetscInt dim1, Nc1, Np1, order1, qa, d1; 2590 PetscInt dim2, Nc2, Np2, order2, qb, d2; 2591 PetscInt dim, Nc, Np, order, qc, d; 2592 2593 PetscFunctionBegin; 2594 PetscValidHeaderSpecific(q1, PETSCQUADRATURE_CLASSID, 1); 2595 PetscValidHeaderSpecific(q2, PETSCQUADRATURE_CLASSID, 2); 2596 PetscAssertPointer(q, 3); 2597 2598 PetscCall(PetscQuadratureGetOrder(q1, &order1)); 2599 PetscCall(PetscQuadratureGetOrder(q2, &order2)); 2600 PetscCheck(order1 == order2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Order1 %" PetscInt_FMT " != %" PetscInt_FMT " Order2", order1, order2); 2601 PetscCall(PetscQuadratureGetData(q1, &dim1, &Nc1, &Np1, &x1, &w1)); 2602 PetscCall(PetscQuadratureGetCellType(q1, &ct1)); 2603 PetscCall(PetscQuadratureGetData(q2, &dim2, &Nc2, &Np2, &x2, &w2)); 2604 PetscCall(PetscQuadratureGetCellType(q2, &ct2)); 2605 PetscCheck(Nc1 == Nc2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "NumComp1 %" PetscInt_FMT " != %" PetscInt_FMT " NumComp2", Nc1, Nc2); 2606 2607 switch (ct1) { 2608 case DM_POLYTOPE_POINT: 2609 ct = ct2; 2610 break; 2611 case DM_POLYTOPE_SEGMENT: 2612 switch (ct2) { 2613 case DM_POLYTOPE_POINT: 2614 ct = ct1; 2615 break; 2616 case DM_POLYTOPE_SEGMENT: 2617 ct = DM_POLYTOPE_QUADRILATERAL; 2618 break; 2619 case DM_POLYTOPE_TRIANGLE: 2620 ct = DM_POLYTOPE_TRI_PRISM; 2621 break; 2622 case DM_POLYTOPE_QUADRILATERAL: 2623 ct = DM_POLYTOPE_HEXAHEDRON; 2624 break; 2625 case DM_POLYTOPE_TETRAHEDRON: 2626 ct = DM_POLYTOPE_UNKNOWN; 2627 break; 2628 case DM_POLYTOPE_HEXAHEDRON: 2629 ct = DM_POLYTOPE_UNKNOWN; 2630 break; 2631 default: 2632 ct = DM_POLYTOPE_UNKNOWN; 2633 } 2634 break; 2635 case DM_POLYTOPE_TRIANGLE: 2636 switch (ct2) { 2637 case DM_POLYTOPE_POINT: 2638 ct = ct1; 2639 break; 2640 case DM_POLYTOPE_SEGMENT: 2641 ct = DM_POLYTOPE_TRI_PRISM; 2642 break; 2643 case DM_POLYTOPE_TRIANGLE: 2644 ct = DM_POLYTOPE_UNKNOWN; 2645 break; 2646 case DM_POLYTOPE_QUADRILATERAL: 2647 ct = DM_POLYTOPE_UNKNOWN; 2648 break; 2649 case DM_POLYTOPE_TETRAHEDRON: 2650 ct = DM_POLYTOPE_UNKNOWN; 2651 break; 2652 case DM_POLYTOPE_HEXAHEDRON: 2653 ct = DM_POLYTOPE_UNKNOWN; 2654 break; 2655 default: 2656 ct = DM_POLYTOPE_UNKNOWN; 2657 } 2658 break; 2659 case DM_POLYTOPE_QUADRILATERAL: 2660 switch (ct2) { 2661 case DM_POLYTOPE_POINT: 2662 ct = ct1; 2663 break; 2664 case DM_POLYTOPE_SEGMENT: 2665 ct = DM_POLYTOPE_HEXAHEDRON; 2666 break; 2667 case DM_POLYTOPE_TRIANGLE: 2668 ct = DM_POLYTOPE_UNKNOWN; 2669 break; 2670 case DM_POLYTOPE_QUADRILATERAL: 2671 ct = DM_POLYTOPE_UNKNOWN; 2672 break; 2673 case DM_POLYTOPE_TETRAHEDRON: 2674 ct = DM_POLYTOPE_UNKNOWN; 2675 break; 2676 case DM_POLYTOPE_HEXAHEDRON: 2677 ct = DM_POLYTOPE_UNKNOWN; 2678 break; 2679 default: 2680 ct = DM_POLYTOPE_UNKNOWN; 2681 } 2682 break; 2683 case DM_POLYTOPE_TETRAHEDRON: 2684 switch (ct2) { 2685 case DM_POLYTOPE_POINT: 2686 ct = ct1; 2687 break; 2688 case DM_POLYTOPE_SEGMENT: 2689 ct = DM_POLYTOPE_UNKNOWN; 2690 break; 2691 case DM_POLYTOPE_TRIANGLE: 2692 ct = DM_POLYTOPE_UNKNOWN; 2693 break; 2694 case DM_POLYTOPE_QUADRILATERAL: 2695 ct = DM_POLYTOPE_UNKNOWN; 2696 break; 2697 case DM_POLYTOPE_TETRAHEDRON: 2698 ct = DM_POLYTOPE_UNKNOWN; 2699 break; 2700 case DM_POLYTOPE_HEXAHEDRON: 2701 ct = DM_POLYTOPE_UNKNOWN; 2702 break; 2703 default: 2704 ct = DM_POLYTOPE_UNKNOWN; 2705 } 2706 break; 2707 case DM_POLYTOPE_HEXAHEDRON: 2708 switch (ct2) { 2709 case DM_POLYTOPE_POINT: 2710 ct = ct1; 2711 break; 2712 case DM_POLYTOPE_SEGMENT: 2713 ct = DM_POLYTOPE_UNKNOWN; 2714 break; 2715 case DM_POLYTOPE_TRIANGLE: 2716 ct = DM_POLYTOPE_UNKNOWN; 2717 break; 2718 case DM_POLYTOPE_QUADRILATERAL: 2719 ct = DM_POLYTOPE_UNKNOWN; 2720 break; 2721 case DM_POLYTOPE_TETRAHEDRON: 2722 ct = DM_POLYTOPE_UNKNOWN; 2723 break; 2724 case DM_POLYTOPE_HEXAHEDRON: 2725 ct = DM_POLYTOPE_UNKNOWN; 2726 break; 2727 default: 2728 ct = DM_POLYTOPE_UNKNOWN; 2729 } 2730 break; 2731 default: 2732 ct = DM_POLYTOPE_UNKNOWN; 2733 } 2734 dim = dim1 + dim2; 2735 Nc = Nc1; 2736 Np = Np1 * Np2; 2737 order = order1; 2738 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 2739 PetscCall(PetscQuadratureSetCellType(*q, ct)); 2740 PetscCall(PetscQuadratureSetOrder(*q, order)); 2741 PetscCall(PetscMalloc1(Np * dim, &x)); 2742 PetscCall(PetscMalloc1(Np, &w)); 2743 for (qa = 0, qc = 0; qa < Np1; ++qa) { 2744 for (qb = 0; qb < Np2; ++qb, ++qc) { 2745 for (d1 = 0, d = 0; d1 < dim1; ++d1, ++d) x[qc * dim + d] = x1[qa * dim1 + d1]; 2746 for (d2 = 0; d2 < dim2; ++d2, ++d) x[qc * dim + d] = x2[qb * dim2 + d2]; 2747 w[qc] = w1[qa] * w2[qb]; 2748 } 2749 } 2750 PetscCall(PetscQuadratureSetData(*q, dim, Nc, Np, x, w)); 2751 PetscFunctionReturn(PETSC_SUCCESS); 2752 } 2753 2754 /* Overwrites A. Can only handle full-rank problems with m>=n 2755 A in column-major format 2756 Ainv in row-major format 2757 tau has length m 2758 worksize must be >= max(1,n) 2759 */ 2760 static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m, PetscInt mstride, PetscInt n, PetscReal *A_in, PetscReal *Ainv_out, PetscScalar *tau, PetscInt worksize, PetscScalar *work) 2761 { 2762 PetscBLASInt M, N, K, lda, ldb, ldwork, info; 2763 PetscScalar *A, *Ainv, *R, *Q, Alpha; 2764 2765 PetscFunctionBegin; 2766 #if defined(PETSC_USE_COMPLEX) 2767 { 2768 PetscInt i, j; 2769 PetscCall(PetscMalloc2(m * n, &A, m * n, &Ainv)); 2770 for (j = 0; j < n; j++) { 2771 for (i = 0; i < m; i++) A[i + m * j] = A_in[i + mstride * j]; 2772 } 2773 mstride = m; 2774 } 2775 #else 2776 A = A_in; 2777 Ainv = Ainv_out; 2778 #endif 2779 2780 PetscCall(PetscBLASIntCast(m, &M)); 2781 PetscCall(PetscBLASIntCast(n, &N)); 2782 PetscCall(PetscBLASIntCast(mstride, &lda)); 2783 PetscCall(PetscBLASIntCast(worksize, &ldwork)); 2784 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 2785 PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info)); 2786 PetscCall(PetscFPTrapPop()); 2787 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error"); 2788 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */ 2789 2790 /* Extract an explicit representation of Q */ 2791 Q = Ainv; 2792 PetscCall(PetscArraycpy(Q, A, mstride * n)); 2793 K = N; /* full rank */ 2794 PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info)); 2795 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error"); 2796 2797 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */ 2798 Alpha = 1.0; 2799 ldb = lda; 2800 PetscCallBLAS("BLAStrsm", BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb)); 2801 /* Ainv is Q, overwritten with inverse */ 2802 2803 #if defined(PETSC_USE_COMPLEX) 2804 { 2805 PetscInt i; 2806 for (i = 0; i < m * n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]); 2807 PetscCall(PetscFree2(A, Ainv)); 2808 } 2809 #endif 2810 PetscFunctionReturn(PETSC_SUCCESS); 2811 } 2812 2813 /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */ 2814 static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval, const PetscReal *x, PetscInt ndegree, const PetscInt *degrees, PetscBool Transpose, PetscReal *B) 2815 { 2816 PetscReal *Bv; 2817 PetscInt i, j; 2818 2819 PetscFunctionBegin; 2820 PetscCall(PetscMalloc1((ninterval + 1) * ndegree, &Bv)); 2821 /* Point evaluation of L_p on all the source vertices */ 2822 PetscCall(PetscDTLegendreEval(ninterval + 1, x, ndegree, degrees, Bv, NULL, NULL)); 2823 /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */ 2824 for (i = 0; i < ninterval; i++) { 2825 for (j = 0; j < ndegree; j++) { 2826 if (Transpose) B[i + ninterval * j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j]; 2827 else B[i * ndegree + j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j]; 2828 } 2829 } 2830 PetscCall(PetscFree(Bv)); 2831 PetscFunctionReturn(PETSC_SUCCESS); 2832 } 2833 2834 /*@ 2835 PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals 2836 2837 Not Collective 2838 2839 Input Parameters: 2840 + degree - degree of reconstruction polynomial 2841 . nsource - number of source intervals 2842 . sourcex - sorted coordinates of source cell boundaries (length `nsource`+1) 2843 . ntarget - number of target intervals 2844 - targetx - sorted coordinates of target cell boundaries (length `ntarget`+1) 2845 2846 Output Parameter: 2847 . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s] 2848 2849 Level: advanced 2850 2851 .seealso: `PetscDTLegendreEval()` 2852 @*/ 2853 PetscErrorCode PetscDTReconstructPoly(PetscInt degree, PetscInt nsource, const PetscReal sourcex[], PetscInt ntarget, const PetscReal targetx[], PetscReal R[]) 2854 { 2855 PetscInt i, j, k, *bdegrees, worksize; 2856 PetscReal xmin, xmax, center, hscale, *sourcey, *targety, *Bsource, *Bsinv, *Btarget; 2857 PetscScalar *tau, *work; 2858 2859 PetscFunctionBegin; 2860 PetscAssertPointer(sourcex, 3); 2861 PetscAssertPointer(targetx, 5); 2862 PetscAssertPointer(R, 6); 2863 PetscCheck(degree < nsource, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Reconstruction degree %" PetscInt_FMT " must be less than number of source intervals %" PetscInt_FMT, degree, nsource); 2864 if (PetscDefined(USE_DEBUG)) { 2865 for (i = 0; i < nsource; i++) PetscCheck(sourcex[i] < sourcex[i + 1], PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Source interval %" PetscInt_FMT " has negative orientation (%g,%g)", i, (double)sourcex[i], (double)sourcex[i + 1]); 2866 for (i = 0; i < ntarget; i++) PetscCheck(targetx[i] < targetx[i + 1], PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Target interval %" PetscInt_FMT " has negative orientation (%g,%g)", i, (double)targetx[i], (double)targetx[i + 1]); 2867 } 2868 xmin = PetscMin(sourcex[0], targetx[0]); 2869 xmax = PetscMax(sourcex[nsource], targetx[ntarget]); 2870 center = (xmin + xmax) / 2; 2871 hscale = (xmax - xmin) / 2; 2872 worksize = nsource; 2873 PetscCall(PetscMalloc4(degree + 1, &bdegrees, nsource + 1, &sourcey, nsource * (degree + 1), &Bsource, worksize, &work)); 2874 PetscCall(PetscMalloc4(nsource, &tau, nsource * (degree + 1), &Bsinv, ntarget + 1, &targety, ntarget * (degree + 1), &Btarget)); 2875 for (i = 0; i <= nsource; i++) sourcey[i] = (sourcex[i] - center) / hscale; 2876 for (i = 0; i <= degree; i++) bdegrees[i] = i + 1; 2877 PetscCall(PetscDTLegendreIntegrate(nsource, sourcey, degree + 1, bdegrees, PETSC_TRUE, Bsource)); 2878 PetscCall(PetscDTPseudoInverseQR(nsource, nsource, degree + 1, Bsource, Bsinv, tau, nsource, work)); 2879 for (i = 0; i <= ntarget; i++) targety[i] = (targetx[i] - center) / hscale; 2880 PetscCall(PetscDTLegendreIntegrate(ntarget, targety, degree + 1, bdegrees, PETSC_FALSE, Btarget)); 2881 for (i = 0; i < ntarget; i++) { 2882 PetscReal rowsum = 0; 2883 for (j = 0; j < nsource; j++) { 2884 PetscReal sum = 0; 2885 for (k = 0; k < degree + 1; k++) sum += Btarget[i * (degree + 1) + k] * Bsinv[k * nsource + j]; 2886 R[i * nsource + j] = sum; 2887 rowsum += sum; 2888 } 2889 for (j = 0; j < nsource; j++) R[i * nsource + j] /= rowsum; /* normalize each row */ 2890 } 2891 PetscCall(PetscFree4(bdegrees, sourcey, Bsource, work)); 2892 PetscCall(PetscFree4(tau, Bsinv, targety, Btarget)); 2893 PetscFunctionReturn(PETSC_SUCCESS); 2894 } 2895 2896 /*@ 2897 PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points 2898 2899 Not Collective 2900 2901 Input Parameters: 2902 + n - the number of GLL nodes 2903 . nodes - the GLL nodes 2904 . weights - the GLL weights 2905 - f - the function values at the nodes 2906 2907 Output Parameter: 2908 . in - the value of the integral 2909 2910 Level: beginner 2911 2912 .seealso: `PetscDTGaussLobattoLegendreQuadrature()` 2913 @*/ 2914 PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n, PetscReal nodes[], PetscReal weights[], const PetscReal f[], PetscReal *in) 2915 { 2916 PetscInt i; 2917 2918 PetscFunctionBegin; 2919 *in = 0.; 2920 for (i = 0; i < n; i++) *in += f[i] * f[i] * weights[i]; 2921 PetscFunctionReturn(PETSC_SUCCESS); 2922 } 2923 2924 /*@C 2925 PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element 2926 2927 Not Collective 2928 2929 Input Parameters: 2930 + n - the number of GLL nodes 2931 . nodes - the GLL nodes, of length `n` 2932 - weights - the GLL weights, of length `n` 2933 2934 Output Parameter: 2935 . AA - the stiffness element, of size `n` by `n` 2936 2937 Level: beginner 2938 2939 Notes: 2940 Destroy this with `PetscGaussLobattoLegendreElementLaplacianDestroy()` 2941 2942 You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row-oriented (the array is symmetric) 2943 2944 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()` 2945 @*/ 2946 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 2947 { 2948 PetscReal **A; 2949 const PetscReal *gllnodes = nodes; 2950 const PetscInt p = n - 1; 2951 PetscReal z0, z1, z2 = -1, x, Lpj, Lpr; 2952 PetscInt i, j, nn, r; 2953 2954 PetscFunctionBegin; 2955 PetscCall(PetscMalloc1(n, &A)); 2956 PetscCall(PetscMalloc1(n * n, &A[0])); 2957 for (i = 1; i < n; i++) A[i] = A[i - 1] + n; 2958 2959 for (j = 1; j < p; j++) { 2960 x = gllnodes[j]; 2961 z0 = 1.; 2962 z1 = x; 2963 for (nn = 1; nn < p; nn++) { 2964 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 2965 z0 = z1; 2966 z1 = z2; 2967 } 2968 Lpj = z2; 2969 for (r = 1; r < p; r++) { 2970 if (r == j) { 2971 A[j][j] = 2. / (3. * (1. - gllnodes[j] * gllnodes[j]) * Lpj * Lpj); 2972 } else { 2973 x = gllnodes[r]; 2974 z0 = 1.; 2975 z1 = x; 2976 for (nn = 1; nn < p; nn++) { 2977 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 2978 z0 = z1; 2979 z1 = z2; 2980 } 2981 Lpr = z2; 2982 A[r][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * Lpr * (gllnodes[j] - gllnodes[r]) * (gllnodes[j] - gllnodes[r])); 2983 } 2984 } 2985 } 2986 for (j = 1; j < p + 1; j++) { 2987 x = gllnodes[j]; 2988 z0 = 1.; 2989 z1 = x; 2990 for (nn = 1; nn < p; nn++) { 2991 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 2992 z0 = z1; 2993 z1 = z2; 2994 } 2995 Lpj = z2; 2996 A[j][0] = 4. * PetscPowRealInt(-1., p) / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. + gllnodes[j]) * (1. + gllnodes[j])); 2997 A[0][j] = A[j][0]; 2998 } 2999 for (j = 0; j < p; j++) { 3000 x = gllnodes[j]; 3001 z0 = 1.; 3002 z1 = x; 3003 for (nn = 1; nn < p; nn++) { 3004 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 3005 z0 = z1; 3006 z1 = z2; 3007 } 3008 Lpj = z2; 3009 3010 A[p][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. - gllnodes[j]) * (1. - gllnodes[j])); 3011 A[j][p] = A[p][j]; 3012 } 3013 A[0][0] = 0.5 + (((PetscReal)p) * (((PetscReal)p) + 1.) - 2.) / 6.; 3014 A[p][p] = A[0][0]; 3015 *AA = A; 3016 PetscFunctionReturn(PETSC_SUCCESS); 3017 } 3018 3019 /*@C 3020 PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element created with `PetscGaussLobattoLegendreElementLaplacianCreate()` 3021 3022 Not Collective 3023 3024 Input Parameters: 3025 + n - the number of GLL nodes 3026 . nodes - the GLL nodes, ignored 3027 . weights - the GLL weightss, ignored 3028 - AA - the stiffness element from `PetscGaussLobattoLegendreElementLaplacianCreate()` 3029 3030 Level: beginner 3031 3032 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()` 3033 @*/ 3034 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 3035 { 3036 PetscFunctionBegin; 3037 PetscCall(PetscFree((*AA)[0])); 3038 PetscCall(PetscFree(*AA)); 3039 *AA = NULL; 3040 PetscFunctionReturn(PETSC_SUCCESS); 3041 } 3042 3043 /*@C 3044 PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element 3045 3046 Not Collective 3047 3048 Input Parameters: 3049 + n - the number of GLL nodes 3050 . nodes - the GLL nodes, of length `n` 3051 - weights - the GLL weights, of length `n` 3052 3053 Output Parameters: 3054 + AA - the stiffness element, of dimension `n` by `n` 3055 - AAT - the transpose of AA (pass in `NULL` if you do not need this array), of dimension `n` by `n` 3056 3057 Level: beginner 3058 3059 Notes: 3060 Destroy this with `PetscGaussLobattoLegendreElementGradientDestroy()` 3061 3062 You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row-oriented 3063 3064 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`, `PetscGaussLobattoLegendreElementGradientDestroy()` 3065 @*/ 3066 PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA, PetscReal ***AAT) 3067 { 3068 PetscReal **A, **AT = NULL; 3069 const PetscReal *gllnodes = nodes; 3070 const PetscInt p = n - 1; 3071 PetscReal Li, Lj, d0; 3072 PetscInt i, j; 3073 3074 PetscFunctionBegin; 3075 PetscCall(PetscMalloc1(n, &A)); 3076 PetscCall(PetscMalloc1(n * n, &A[0])); 3077 for (i = 1; i < n; i++) A[i] = A[i - 1] + n; 3078 3079 if (AAT) { 3080 PetscCall(PetscMalloc1(n, &AT)); 3081 PetscCall(PetscMalloc1(n * n, &AT[0])); 3082 for (i = 1; i < n; i++) AT[i] = AT[i - 1] + n; 3083 } 3084 3085 if (n == 1) A[0][0] = 0.; 3086 d0 = (PetscReal)p * ((PetscReal)p + 1.) / 4.; 3087 for (i = 0; i < n; i++) { 3088 for (j = 0; j < n; j++) { 3089 A[i][j] = 0.; 3090 PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li)); 3091 PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj)); 3092 if (i != j) A[i][j] = Li / (Lj * (gllnodes[i] - gllnodes[j])); 3093 if ((j == i) && (i == 0)) A[i][j] = -d0; 3094 if (j == i && i == p) A[i][j] = d0; 3095 if (AT) AT[j][i] = A[i][j]; 3096 } 3097 } 3098 if (AAT) *AAT = AT; 3099 *AA = A; 3100 PetscFunctionReturn(PETSC_SUCCESS); 3101 } 3102 3103 /*@C 3104 PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with `PetscGaussLobattoLegendreElementGradientCreate()` 3105 3106 Not Collective 3107 3108 Input Parameters: 3109 + n - the number of GLL nodes 3110 . nodes - the GLL nodes, ignored 3111 . weights - the GLL weights, ignored 3112 . AA - the stiffness element obtained with `PetscGaussLobattoLegendreElementGradientCreate()` 3113 - AAT - the transpose of the element obtained with `PetscGaussLobattoLegendreElementGradientCreate()` 3114 3115 Level: beginner 3116 3117 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionCreate()` 3118 @*/ 3119 PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA, PetscReal ***AAT) 3120 { 3121 PetscFunctionBegin; 3122 PetscCall(PetscFree((*AA)[0])); 3123 PetscCall(PetscFree(*AA)); 3124 *AA = NULL; 3125 if (AAT) { 3126 PetscCall(PetscFree((*AAT)[0])); 3127 PetscCall(PetscFree(*AAT)); 3128 *AAT = NULL; 3129 } 3130 PetscFunctionReturn(PETSC_SUCCESS); 3131 } 3132 3133 /*@C 3134 PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element 3135 3136 Not Collective 3137 3138 Input Parameters: 3139 + n - the number of GLL nodes 3140 . nodes - the GLL nodes, of length `n` 3141 - weights - the GLL weights, of length `n` 3142 3143 Output Parameter: 3144 . AA - the stiffness element, of dimension `n` by `n` 3145 3146 Level: beginner 3147 3148 Notes: 3149 Destroy this with `PetscGaussLobattoLegendreElementAdvectionDestroy()` 3150 3151 This is the same as the Gradient operator multiplied by the diagonal mass matrix 3152 3153 You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row-oriented 3154 3155 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionDestroy()` 3156 @*/ 3157 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 3158 { 3159 PetscReal **D; 3160 const PetscReal *gllweights = weights; 3161 const PetscInt glln = n; 3162 PetscInt i, j; 3163 3164 PetscFunctionBegin; 3165 PetscCall(PetscGaussLobattoLegendreElementGradientCreate(n, nodes, weights, &D, NULL)); 3166 for (i = 0; i < glln; i++) { 3167 for (j = 0; j < glln; j++) D[i][j] = gllweights[i] * D[i][j]; 3168 } 3169 *AA = D; 3170 PetscFunctionReturn(PETSC_SUCCESS); 3171 } 3172 3173 /*@C 3174 PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element created with `PetscGaussLobattoLegendreElementAdvectionCreate()` 3175 3176 Not Collective 3177 3178 Input Parameters: 3179 + n - the number of GLL nodes 3180 . nodes - the GLL nodes, ignored 3181 . weights - the GLL weights, ignored 3182 - AA - advection obtained with `PetscGaussLobattoLegendreElementAdvectionCreate()` 3183 3184 Level: beginner 3185 3186 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementAdvectionCreate()` 3187 @*/ 3188 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 3189 { 3190 PetscFunctionBegin; 3191 PetscCall(PetscFree((*AA)[0])); 3192 PetscCall(PetscFree(*AA)); 3193 *AA = NULL; 3194 PetscFunctionReturn(PETSC_SUCCESS); 3195 } 3196 3197 PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 3198 { 3199 PetscReal **A; 3200 const PetscReal *gllweights = weights; 3201 const PetscInt glln = n; 3202 PetscInt i, j; 3203 3204 PetscFunctionBegin; 3205 PetscCall(PetscMalloc1(glln, &A)); 3206 PetscCall(PetscMalloc1(glln * glln, &A[0])); 3207 for (i = 1; i < glln; i++) A[i] = A[i - 1] + glln; 3208 if (glln == 1) A[0][0] = 0.; 3209 for (i = 0; i < glln; i++) { 3210 for (j = 0; j < glln; j++) { 3211 A[i][j] = 0.; 3212 if (j == i) A[i][j] = gllweights[i]; 3213 } 3214 } 3215 *AA = A; 3216 PetscFunctionReturn(PETSC_SUCCESS); 3217 } 3218 3219 PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 3220 { 3221 PetscFunctionBegin; 3222 PetscCall(PetscFree((*AA)[0])); 3223 PetscCall(PetscFree(*AA)); 3224 *AA = NULL; 3225 PetscFunctionReturn(PETSC_SUCCESS); 3226 } 3227 3228 /*@ 3229 PetscDTIndexToBary - convert an index into a barycentric coordinate. 3230 3231 Input Parameters: 3232 + len - the desired length of the barycentric tuple (usually 1 more than the dimension it represents, so a barycentric coordinate in a triangle has length 3) 3233 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 3234 - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum) 3235 3236 Output Parameter: 3237 . coord - will be filled with the barycentric coordinate, of length `n` 3238 3239 Level: beginner 3240 3241 Note: 3242 The indices map to barycentric coordinates in lexicographic order, where the first index is the 3243 least significant and the last index is the most significant. 3244 3245 .seealso: `PetscDTBaryToIndex()` 3246 @*/ 3247 PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[]) 3248 { 3249 PetscInt c, d, s, total, subtotal, nexttotal; 3250 3251 PetscFunctionBeginHot; 3252 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 3253 PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative"); 3254 if (!len) { 3255 if (!sum && !index) PetscFunctionReturn(PETSC_SUCCESS); 3256 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 3257 } 3258 for (c = 1, total = 1; c <= len; c++) { 3259 /* total is the number of ways to have a tuple of length c with sum */ 3260 if (index < total) break; 3261 total = (total * (sum + c)) / c; 3262 } 3263 PetscCheck(c <= len, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range"); 3264 for (d = c; d < len; d++) coord[d] = 0; 3265 for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) { 3266 /* subtotal is the number of ways to have a tuple of length c with sum s */ 3267 /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */ 3268 if ((index + subtotal) >= total) { 3269 coord[--c] = sum - s; 3270 index -= (total - subtotal); 3271 sum = s; 3272 total = nexttotal; 3273 subtotal = 1; 3274 nexttotal = 1; 3275 s = 0; 3276 } else { 3277 subtotal = (subtotal * (c + s)) / (s + 1); 3278 nexttotal = (nexttotal * (c - 1 + s)) / (s + 1); 3279 s++; 3280 } 3281 } 3282 PetscFunctionReturn(PETSC_SUCCESS); 3283 } 3284 3285 /*@ 3286 PetscDTBaryToIndex - convert a barycentric coordinate to an index 3287 3288 Input Parameters: 3289 + len - the desired length of the barycentric tuple (usually 1 more than the dimension it represents, so a barycentric coordinate in a triangle has length 3) 3290 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 3291 - coord - a barycentric coordinate with the given length `len` and `sum` 3292 3293 Output Parameter: 3294 . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum) 3295 3296 Level: beginner 3297 3298 Note: 3299 The indices map to barycentric coordinates in lexicographic order, where the first index is the 3300 least significant and the last index is the most significant. 3301 3302 .seealso: `PetscDTIndexToBary` 3303 @*/ 3304 PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index) 3305 { 3306 PetscInt c; 3307 PetscInt i; 3308 PetscInt total; 3309 3310 PetscFunctionBeginHot; 3311 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 3312 if (!len) { 3313 if (!sum) { 3314 *index = 0; 3315 PetscFunctionReturn(PETSC_SUCCESS); 3316 } 3317 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 3318 } 3319 for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c; 3320 i = total - 1; 3321 c = len - 1; 3322 sum -= coord[c]; 3323 while (sum > 0) { 3324 PetscInt subtotal; 3325 PetscInt s; 3326 3327 for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s; 3328 i -= subtotal; 3329 sum -= coord[--c]; 3330 } 3331 *index = i; 3332 PetscFunctionReturn(PETSC_SUCCESS); 3333 } 3334 3335 /*@ 3336 PetscQuadratureComputePermutations - Compute permutations of quadrature points corresponding to domain orientations 3337 3338 Input Parameter: 3339 . quad - The `PetscQuadrature` 3340 3341 Output Parameters: 3342 + Np - The number of domain orientations 3343 - perm - An array of `IS` permutations, one for ech orientation, 3344 3345 Level: developer 3346 3347 .seealso: `PetscQuadratureSetCellType()`, `PetscQuadrature` 3348 @*/ 3349 PetscErrorCode PetscQuadratureComputePermutations(PetscQuadrature quad, PeOp PetscInt *Np, IS *perm[]) 3350 { 3351 DMPolytopeType ct; 3352 const PetscReal *xq, *wq; 3353 PetscInt dim, qdim, d, Na, o, Nq, q, qp; 3354 3355 PetscFunctionBegin; 3356 PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &xq, &wq)); 3357 PetscCall(PetscQuadratureGetCellType(quad, &ct)); 3358 dim = DMPolytopeTypeGetDim(ct); 3359 Na = DMPolytopeTypeGetNumArrangements(ct); 3360 PetscCall(PetscMalloc1(Na, perm)); 3361 if (Np) *Np = Na; 3362 Na /= 2; 3363 for (o = -Na; o < Na; ++o) { 3364 DM refdm; 3365 PetscInt *idx; 3366 PetscReal xi0[3] = {-1., -1., -1.}, v0[3], J[9], detJ, txq[3]; 3367 PetscBool flg; 3368 3369 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm)); 3370 PetscCall(DMPlexOrientPoint(refdm, 0, o)); 3371 PetscCall(DMPlexComputeCellGeometryFEM(refdm, 0, NULL, v0, J, NULL, &detJ)); 3372 PetscCall(PetscMalloc1(Nq, &idx)); 3373 for (q = 0; q < Nq; ++q) { 3374 CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q * dim], txq); 3375 for (qp = 0; qp < Nq; ++qp) { 3376 PetscReal diff = 0.; 3377 3378 for (d = 0; d < dim; ++d) diff += PetscAbsReal(txq[d] - xq[qp * dim + d]); 3379 if (diff < PETSC_SMALL) break; 3380 } 3381 PetscCheck(qp < Nq, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Transformed quad point %" PetscInt_FMT " does not match another quad point", q); 3382 idx[q] = qp; 3383 } 3384 PetscCall(DMDestroy(&refdm)); 3385 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nq, idx, PETSC_OWN_POINTER, &(*perm)[o + Na])); 3386 PetscCall(ISGetInfo((*perm)[o + Na], IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg)); 3387 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ordering for orientation %" PetscInt_FMT " was not a permutation", o); 3388 PetscCall(ISSetPermutation((*perm)[o + Na])); 3389 } 3390 if (!Na) (*perm)[0] = NULL; 3391 PetscFunctionReturn(PETSC_SUCCESS); 3392 } 3393 3394 /*@ 3395 PetscDTCreateDefaultQuadrature - Create default quadrature for a given cell 3396 3397 Not collective 3398 3399 Input Parameters: 3400 + ct - The integration domain 3401 - qorder - The desired quadrature order 3402 3403 Output Parameters: 3404 + q - The cell quadrature 3405 - fq - The face quadrature 3406 3407 Level: developer 3408 3409 .seealso: `PetscFECreateDefault()`, `PetscDTGaussTensorQuadrature()`, `PetscDTSimplexQuadrature()`, `PetscDTTensorQuadratureCreate()` 3410 @*/ 3411 PetscErrorCode PetscDTCreateDefaultQuadrature(DMPolytopeType ct, PetscInt qorder, PetscQuadrature *q, PetscQuadrature *fq) 3412 { 3413 const PetscInt quadPointsPerEdge = PetscMax(qorder + 1, 1); 3414 const PetscInt dim = DMPolytopeTypeGetDim(ct); 3415 3416 PetscFunctionBegin; 3417 switch (ct) { 3418 case DM_POLYTOPE_SEGMENT: 3419 case DM_POLYTOPE_POINT_PRISM_TENSOR: 3420 case DM_POLYTOPE_QUADRILATERAL: 3421 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3422 case DM_POLYTOPE_HEXAHEDRON: 3423 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3424 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q)); 3425 PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq)); 3426 break; 3427 case DM_POLYTOPE_TRIANGLE: 3428 case DM_POLYTOPE_TETRAHEDRON: 3429 PetscCall(PetscDTSimplexQuadrature(dim, 2 * qorder, PETSCDTSIMPLEXQUAD_DEFAULT, q)); 3430 PetscCall(PetscDTSimplexQuadrature(dim - 1, 2 * qorder, PETSCDTSIMPLEXQUAD_DEFAULT, fq)); 3431 break; 3432 case DM_POLYTOPE_TRI_PRISM: 3433 case DM_POLYTOPE_TRI_PRISM_TENSOR: { 3434 PetscQuadrature q1, q2; 3435 3436 // TODO: this should be able to use symmetric rules, but doing so causes tests to fail 3437 PetscCall(PetscDTSimplexQuadrature(2, 2 * qorder, PETSCDTSIMPLEXQUAD_CONIC, &q1)); 3438 PetscCall(PetscDTGaussTensorQuadrature(1, 1, quadPointsPerEdge, -1.0, 1.0, &q2)); 3439 PetscCall(PetscDTTensorQuadratureCreate(q1, q2, q)); 3440 PetscCall(PetscQuadratureDestroy(&q2)); 3441 *fq = q1; 3442 /* TODO Need separate quadratures for each face */ 3443 } break; 3444 default: 3445 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No quadrature for celltype %s", DMPolytopeTypes[PetscMin(ct, DM_POLYTOPE_UNKNOWN)]); 3446 } 3447 PetscFunctionReturn(PETSC_SUCCESS); 3448 } 3449