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", "PetscDTNodeType", "PETSCDTNODES_", NULL}; 17 const char *const *const PetscDTNodeTypes = PetscDTNodeTypes_shifted + 1; 18 19 const char *const PetscDTSimplexQuadratureTypes_shifted[] = {"default", "conic", "minsym", "diagsym", "PetscDTSimplexQuadratureType", "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 static PetscBool DiagSymTriQuadCite = PETSC_FALSE; 2096 const char DiagSymTriQuadCitation[] = "@article{KongMulderVeldhuizen1999,\n" 2097 " title = {Higher-order triangular and tetrahedral finite elements with mass lumping for solving the wave equation},\n" 2098 " journal = {Journal of Engineering Mathematics},\n" 2099 " volume = {35},\n" 2100 " number = {4},\n" 2101 " pages = {405--426},\n" 2102 " year = {1999},\n" 2103 " doi = {10.1023/A:1004420829610},\n" 2104 " url = {https://link.springer.com/article/10.1023/A:1004420829610},\n" 2105 " author = {MJS Chin-Joe-Kong and Wim A Mulder and Marinus Van Veldhuizen},\n" 2106 "}\n"; 2107 2108 #include "petscdttridiagquadrules.h" 2109 2110 // https://en.wikipedia.org/wiki/Partition_(number_theory) 2111 static PetscErrorCode PetscDTPartitionNumber(PetscInt n, PetscInt *p) 2112 { 2113 // sequence A000041 in the OEIS 2114 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}; 2115 PetscInt tabulated_max = PETSC_STATIC_ARRAY_LENGTH(partition) - 1; 2116 2117 PetscFunctionBegin; 2118 PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Partition number not defined for negative number %" PetscInt_FMT, n); 2119 // not implementing the pentagonal number recurrence, we don't need partition numbers for n that high 2120 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); 2121 *p = partition[n]; 2122 PetscFunctionReturn(PETSC_SUCCESS); 2123 } 2124 2125 /*@ 2126 PetscDTSimplexQuadrature - Create a quadrature rule for a simplex that exactly integrates polynomials up to a given degree. 2127 2128 Not Collective 2129 2130 Input Parameters: 2131 + dim - The spatial dimension of the simplex (1 = segment, 2 = triangle, 3 = tetrahedron) 2132 . degree - The largest polynomial degree that is required to be integrated exactly 2133 - type - `PetscDTSimplexQuadratureType` indicating the type of quadrature rule 2134 2135 Output Parameter: 2136 . quad - A `PetscQuadrature` object for integration over the biunit simplex 2137 (defined by the bounds $x_i >= -1$ and $\sum_i x_i <= 2 - d$) that is exact for 2138 polynomials up to the given degree 2139 2140 Level: intermediate 2141 2142 .seealso: `PetscDTSimplexQuadratureType`, `PetscDTGaussQuadrature()`, `PetscDTStroudCononicalQuadrature()`, `PetscQuadrature` 2143 @*/ 2144 PetscErrorCode PetscDTSimplexQuadrature(PetscInt dim, PetscInt degree, PetscDTSimplexQuadratureType type, PetscQuadrature *quad) 2145 { 2146 PetscDTSimplexQuadratureType orig_type = type; 2147 2148 PetscFunctionBegin; 2149 PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative dimension %" PetscInt_FMT, dim); 2150 PetscCheck(degree >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %" PetscInt_FMT, degree); 2151 if (type == PETSCDTSIMPLEXQUAD_DEFAULT) type = PETSCDTSIMPLEXQUAD_MINSYM; 2152 if (type == PETSCDTSIMPLEXQUAD_CONIC || dim < 2) { 2153 PetscInt points_per_dim = (degree + 2) / 2; // ceil((degree + 1) / 2); 2154 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, points_per_dim, -1, 1, quad)); 2155 } else { 2156 DMPolytopeType ct; 2157 PetscInt n = dim + 1; 2158 PetscInt fact = 1; 2159 PetscInt *part, *perm; 2160 PetscInt p = 0; 2161 PetscInt max_degree; 2162 const PetscInt *nodes_per_type = NULL; 2163 const PetscInt *all_num_full_nodes = NULL; 2164 const PetscReal **weights_list = NULL; 2165 const PetscReal **compact_nodes_list = NULL; 2166 const char *citation = NULL; 2167 PetscBool *cited = NULL; 2168 2169 switch (dim) { 2170 case 0: 2171 ct = DM_POLYTOPE_POINT; 2172 break; 2173 case 1: 2174 ct = DM_POLYTOPE_SEGMENT; 2175 break; 2176 case 2: 2177 ct = DM_POLYTOPE_TRIANGLE; 2178 break; 2179 case 3: 2180 ct = DM_POLYTOPE_TETRAHEDRON; 2181 break; 2182 default: 2183 ct = DM_POLYTOPE_UNKNOWN; 2184 } 2185 if (type == PETSCDTSIMPLEXQUAD_MINSYM) { 2186 switch (dim) { 2187 case 2: 2188 cited = &MinSymTriQuadCite; 2189 citation = MinSymTriQuadCitation; 2190 max_degree = PetscDTWVTriQuad_max_degree; 2191 nodes_per_type = PetscDTWVTriQuad_num_orbits; 2192 all_num_full_nodes = PetscDTWVTriQuad_num_nodes; 2193 weights_list = PetscDTWVTriQuad_weights; 2194 compact_nodes_list = PetscDTWVTriQuad_orbits; 2195 break; 2196 case 3: 2197 cited = &MinSymTetQuadCite; 2198 citation = MinSymTetQuadCitation; 2199 max_degree = PetscDTJSTetQuad_max_degree; 2200 nodes_per_type = PetscDTJSTetQuad_num_orbits; 2201 all_num_full_nodes = PetscDTJSTetQuad_num_nodes; 2202 weights_list = PetscDTJSTetQuad_weights; 2203 compact_nodes_list = PetscDTJSTetQuad_orbits; 2204 break; 2205 default: 2206 max_degree = -1; 2207 break; 2208 } 2209 } else { 2210 switch (dim) { 2211 case 2: 2212 cited = &DiagSymTriQuadCite; 2213 citation = DiagSymTriQuadCitation; 2214 max_degree = PetscDTKMVTriQuad_max_degree; 2215 nodes_per_type = PetscDTKMVTriQuad_num_orbits; 2216 all_num_full_nodes = PetscDTKMVTriQuad_num_nodes; 2217 weights_list = PetscDTKMVTriQuad_weights; 2218 compact_nodes_list = PetscDTKMVTriQuad_orbits; 2219 break; 2220 default: 2221 max_degree = -1; 2222 break; 2223 } 2224 } 2225 2226 if (degree > max_degree) { 2227 if (orig_type == PETSCDTSIMPLEXQUAD_DEFAULT) { 2228 // fall back to conic 2229 PetscCall(PetscDTSimplexQuadrature(dim, degree, PETSCDTSIMPLEXQUAD_CONIC, quad)); 2230 PetscFunctionReturn(PETSC_SUCCESS); 2231 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "%s symmetric quadrature for dim %" PetscInt_FMT ", degree %" PetscInt_FMT " unsupported", orig_type == PETSCDTSIMPLEXQUAD_MINSYM ? "Minimal" : "Diagonal", dim, degree); 2232 } 2233 2234 PetscCall(PetscCitationsRegister(citation, cited)); 2235 2236 PetscCall(PetscDTPartitionNumber(n, &p)); 2237 for (PetscInt d = 2; d <= n; d++) fact *= d; 2238 2239 PetscInt num_full_nodes = all_num_full_nodes[degree]; 2240 const PetscReal *all_compact_nodes = compact_nodes_list[degree]; 2241 const PetscReal *all_compact_weights = weights_list[degree]; 2242 nodes_per_type = &nodes_per_type[p * degree]; 2243 2244 PetscReal *points; 2245 PetscReal *counts; 2246 PetscReal *weights; 2247 PetscReal *bary_to_biunit; // row-major transformation of barycentric coordinate to biunit 2248 PetscQuadrature q; 2249 2250 // compute the transformation 2251 PetscCall(PetscMalloc1(n * dim, &bary_to_biunit)); 2252 for (PetscInt d = 0; d < dim; d++) { 2253 for (PetscInt b = 0; b < n; b++) bary_to_biunit[d * n + b] = (d == b) ? 1.0 : -1.0; 2254 } 2255 2256 PetscCall(PetscMalloc3(n, &part, n, &perm, n, &counts)); 2257 PetscCall(PetscCalloc1(num_full_nodes * dim, &points)); 2258 PetscCall(PetscMalloc1(num_full_nodes, &weights)); 2259 2260 // (0, 0, ...) is the first partition lexicographically 2261 PetscCall(PetscArrayzero(part, n)); 2262 PetscCall(PetscArrayzero(counts, n)); 2263 counts[0] = n; 2264 2265 // for each partition 2266 for (PetscInt s = 0, node_offset = 0; s < p; s++) { 2267 PetscInt num_compact_coords = part[n - 1] + 1; 2268 2269 const PetscReal *compact_nodes = all_compact_nodes; 2270 const PetscReal *compact_weights = all_compact_weights; 2271 all_compact_nodes += num_compact_coords * nodes_per_type[s]; 2272 all_compact_weights += nodes_per_type[s]; 2273 2274 // for every permutation of the vertices 2275 for (PetscInt f = 0; f < fact; f++) { 2276 PetscCall(PetscDTEnumPerm(n, f, perm, NULL)); 2277 2278 // check if it is a valid permutation 2279 PetscInt digit; 2280 for (digit = 1; digit < n; digit++) { 2281 // skip permutations that would duplicate a node because it has a smaller symmetry group 2282 if (part[digit - 1] == part[digit] && perm[digit - 1] > perm[digit]) break; 2283 } 2284 if (digit < n) continue; 2285 2286 // create full nodes from this permutation of the compact nodes 2287 PetscReal *full_nodes = &points[node_offset * dim]; 2288 PetscReal *full_weights = &weights[node_offset]; 2289 2290 PetscCall(PetscArraycpy(full_weights, compact_weights, nodes_per_type[s])); 2291 for (PetscInt b = 0; b < n; b++) { 2292 for (PetscInt d = 0; d < dim; d++) { 2293 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]]; 2294 } 2295 } 2296 node_offset += nodes_per_type[s]; 2297 } 2298 2299 if (s < p - 1) { // Generate the next partition 2300 /* A partition is described by the number of coordinates that are in 2301 * each set of duplicates (counts) and redundantly by mapping each 2302 * index to its set of duplicates (part) 2303 * 2304 * Counts should always be in nonincreasing order 2305 * 2306 * We want to generate the partitions lexically by part, which means 2307 * finding the last index where count > 1 and reducing by 1. 2308 * 2309 * For the new counts beyond that index, we eagerly assign the remaining 2310 * capacity of the partition to smaller indices (ensures lexical ordering), 2311 * while respecting the nonincreasing invariant of the counts 2312 */ 2313 PetscInt last_digit = part[n - 1]; 2314 PetscInt last_digit_with_extra = last_digit; 2315 while (counts[last_digit_with_extra] == 1) last_digit_with_extra--; 2316 PetscInt limit = --counts[last_digit_with_extra]; 2317 PetscInt total_to_distribute = last_digit - last_digit_with_extra + 1; 2318 for (PetscInt digit = last_digit_with_extra + 1; digit < n; digit++) { 2319 counts[digit] = PetscMin(limit, total_to_distribute); 2320 total_to_distribute -= PetscMin(limit, total_to_distribute); 2321 } 2322 for (PetscInt digit = 0, offset = 0; digit < n; digit++) { 2323 PetscInt count = counts[digit]; 2324 for (PetscInt c = 0; c < count; c++) part[offset++] = digit; 2325 } 2326 } 2327 PetscCheck(node_offset <= num_full_nodes, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node offset %" PetscInt_FMT " > %" PetscInt_FMT " number of nodes", node_offset, num_full_nodes); 2328 } 2329 PetscCall(PetscFree3(part, perm, counts)); 2330 PetscCall(PetscFree(bary_to_biunit)); 2331 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &q)); 2332 PetscCall(PetscQuadratureSetCellType(q, ct)); 2333 PetscCall(PetscQuadratureSetOrder(q, degree)); 2334 PetscCall(PetscQuadratureSetData(q, dim, 1, num_full_nodes, points, weights)); 2335 *quad = q; 2336 } 2337 PetscFunctionReturn(PETSC_SUCCESS); 2338 } 2339 2340 /*@ 2341 PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell 2342 2343 Not Collective 2344 2345 Input Parameters: 2346 + dim - The cell dimension 2347 . level - The number of points in one dimension, $2^l$ 2348 . a - left end of interval (often-1) 2349 - b - right end of interval (often +1) 2350 2351 Output Parameter: 2352 . q - A `PetscQuadrature` object 2353 2354 Level: intermediate 2355 2356 .seealso: `PetscDTGaussTensorQuadrature()`, `PetscQuadrature` 2357 @*/ 2358 PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q) 2359 { 2360 DMPolytopeType ct; 2361 const PetscInt p = 16; /* Digits of precision in the evaluation */ 2362 const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */ 2363 const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */ 2364 const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */ 2365 PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */ 2366 PetscReal wk = 0.5 * PETSC_PI; /* Quadrature weight at x_k */ 2367 PetscReal *x, *w; 2368 PetscInt K, k, npoints; 2369 2370 PetscFunctionBegin; 2371 PetscCheck(dim <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not yet implemented", dim); 2372 PetscCheck(level, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits"); 2373 switch (dim) { 2374 case 0: 2375 ct = DM_POLYTOPE_POINT; 2376 break; 2377 case 1: 2378 ct = DM_POLYTOPE_SEGMENT; 2379 break; 2380 case 2: 2381 ct = DM_POLYTOPE_QUADRILATERAL; 2382 break; 2383 case 3: 2384 ct = DM_POLYTOPE_HEXAHEDRON; 2385 break; 2386 default: 2387 ct = DM_POLYTOPE_UNKNOWN; 2388 } 2389 /* Find K such that the weights are < 32 digits of precision */ 2390 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))); 2391 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 2392 PetscCall(PetscQuadratureSetCellType(*q, ct)); 2393 PetscCall(PetscQuadratureSetOrder(*q, 2 * K + 1)); 2394 npoints = 2 * K - 1; 2395 PetscCall(PetscMalloc1(npoints * dim, &x)); 2396 PetscCall(PetscMalloc1(npoints, &w)); 2397 /* Center term */ 2398 x[0] = beta; 2399 w[0] = 0.5 * alpha * PETSC_PI; 2400 for (k = 1; k < K; ++k) { 2401 wk = 0.5 * alpha * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2402 xk = PetscTanhReal(0.5 * PETSC_PI * PetscSinhReal(k * h)); 2403 x[2 * k - 1] = -alpha * xk + beta; 2404 w[2 * k - 1] = wk; 2405 x[2 * k + 0] = alpha * xk + beta; 2406 w[2 * k + 0] = wk; 2407 } 2408 PetscCall(PetscQuadratureSetData(*q, dim, 1, npoints, x, w)); 2409 PetscFunctionReturn(PETSC_SUCCESS); 2410 } 2411 2412 PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2413 { 2414 const PetscInt p = 16; /* Digits of precision in the evaluation */ 2415 const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */ 2416 const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */ 2417 PetscReal h = 1.0; /* Step size, length between x_k */ 2418 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 2419 PetscReal osum = 0.0; /* Integral on last level */ 2420 PetscReal psum = 0.0; /* Integral on the level before the last level */ 2421 PetscReal sum; /* Integral on current level */ 2422 PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 2423 PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 2424 PetscReal wk; /* Quadrature weight at x_k */ 2425 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 2426 PetscInt d; /* Digits of precision in the integral */ 2427 2428 PetscFunctionBegin; 2429 PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 2430 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 2431 /* Center term */ 2432 func(&beta, ctx, &lval); 2433 sum = 0.5 * alpha * PETSC_PI * lval; 2434 /* */ 2435 do { 2436 PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4; 2437 PetscInt k = 1; 2438 2439 ++l; 2440 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */ 2441 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 2442 psum = osum; 2443 osum = sum; 2444 h *= 0.5; 2445 sum *= 0.5; 2446 do { 2447 wk = 0.5 * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2448 yk = 1.0 / (PetscExpReal(0.5 * PETSC_PI * PetscSinhReal(k * h)) * PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2449 lx = -alpha * (1.0 - yk) + beta; 2450 rx = alpha * (1.0 - yk) + beta; 2451 func(&lx, ctx, &lval); 2452 func(&rx, ctx, &rval); 2453 lterm = alpha * wk * lval; 2454 maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm); 2455 sum += lterm; 2456 rterm = alpha * wk * rval; 2457 maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm); 2458 sum += rterm; 2459 ++k; 2460 /* Only need to evaluate every other point on refined levels */ 2461 if (l != 1) ++k; 2462 } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */ 2463 2464 d1 = PetscLog10Real(PetscAbsReal(sum - osum)); 2465 d2 = PetscLog10Real(PetscAbsReal(sum - psum)); 2466 d3 = PetscLog10Real(maxTerm) - p; 2467 if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0; 2468 else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm))); 2469 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4))); 2470 } while (d < digits && l < 12); 2471 *sol = sum; 2472 PetscCall(PetscFPTrapPop()); 2473 PetscFunctionReturn(PETSC_SUCCESS); 2474 } 2475 2476 #if defined(PETSC_HAVE_MPFR) 2477 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2478 { 2479 const PetscInt safetyFactor = 2; /* Calculate abscissa until 2*p digits */ 2480 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 2481 mpfr_t alpha; /* Half-width of the integration interval */ 2482 mpfr_t beta; /* Center of the integration interval */ 2483 mpfr_t h; /* Step size, length between x_k */ 2484 mpfr_t osum; /* Integral on last level */ 2485 mpfr_t psum; /* Integral on the level before the last level */ 2486 mpfr_t sum; /* Integral on current level */ 2487 mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 2488 mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 2489 mpfr_t wk; /* Quadrature weight at x_k */ 2490 PetscReal lval, rval, rtmp; /* Terms in the quadature sum to the left and right of 0 */ 2491 PetscInt d; /* Digits of precision in the integral */ 2492 mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp; 2493 2494 PetscFunctionBegin; 2495 PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 2496 /* Create high precision storage */ 2497 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); 2498 /* Initialization */ 2499 mpfr_set_d(alpha, 0.5 * (b - a), MPFR_RNDN); 2500 mpfr_set_d(beta, 0.5 * (b + a), MPFR_RNDN); 2501 mpfr_set_d(osum, 0.0, MPFR_RNDN); 2502 mpfr_set_d(psum, 0.0, MPFR_RNDN); 2503 mpfr_set_d(h, 1.0, MPFR_RNDN); 2504 mpfr_const_pi(pi2, MPFR_RNDN); 2505 mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN); 2506 /* Center term */ 2507 rtmp = 0.5 * (b + a); 2508 func(&rtmp, ctx, &lval); 2509 mpfr_set(sum, pi2, MPFR_RNDN); 2510 mpfr_mul(sum, sum, alpha, MPFR_RNDN); 2511 mpfr_mul_d(sum, sum, lval, MPFR_RNDN); 2512 /* */ 2513 do { 2514 PetscReal d1, d2, d3, d4; 2515 PetscInt k = 1; 2516 2517 ++l; 2518 mpfr_set_d(maxTerm, 0.0, MPFR_RNDN); 2519 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */ 2520 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 2521 mpfr_set(psum, osum, MPFR_RNDN); 2522 mpfr_set(osum, sum, MPFR_RNDN); 2523 mpfr_mul_d(h, h, 0.5, MPFR_RNDN); 2524 mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN); 2525 do { 2526 mpfr_set_si(kh, k, MPFR_RNDN); 2527 mpfr_mul(kh, kh, h, MPFR_RNDN); 2528 /* Weight */ 2529 mpfr_set(wk, h, MPFR_RNDN); 2530 mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN); 2531 mpfr_mul(msinh, msinh, pi2, MPFR_RNDN); 2532 mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN); 2533 mpfr_cosh(tmp, msinh, MPFR_RNDN); 2534 mpfr_sqr(tmp, tmp, MPFR_RNDN); 2535 mpfr_mul(wk, wk, mcosh, MPFR_RNDN); 2536 mpfr_div(wk, wk, tmp, MPFR_RNDN); 2537 /* Abscissa */ 2538 mpfr_set_d(yk, 1.0, MPFR_RNDZ); 2539 mpfr_cosh(tmp, msinh, MPFR_RNDN); 2540 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 2541 mpfr_exp(tmp, msinh, MPFR_RNDN); 2542 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 2543 /* Quadrature points */ 2544 mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ); 2545 mpfr_mul(lx, lx, alpha, MPFR_RNDU); 2546 mpfr_add(lx, lx, beta, MPFR_RNDU); 2547 mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ); 2548 mpfr_mul(rx, rx, alpha, MPFR_RNDD); 2549 mpfr_add(rx, rx, beta, MPFR_RNDD); 2550 /* Evaluation */ 2551 rtmp = mpfr_get_d(lx, MPFR_RNDU); 2552 func(&rtmp, ctx, &lval); 2553 rtmp = mpfr_get_d(rx, MPFR_RNDD); 2554 func(&rtmp, ctx, &rval); 2555 /* Update */ 2556 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 2557 mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN); 2558 mpfr_add(sum, sum, tmp, MPFR_RNDN); 2559 mpfr_abs(tmp, tmp, MPFR_RNDN); 2560 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 2561 mpfr_set(curTerm, tmp, MPFR_RNDN); 2562 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 2563 mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN); 2564 mpfr_add(sum, sum, tmp, MPFR_RNDN); 2565 mpfr_abs(tmp, tmp, MPFR_RNDN); 2566 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 2567 mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN); 2568 ++k; 2569 /* Only need to evaluate every other point on refined levels */ 2570 if (l != 1) ++k; 2571 mpfr_log10(tmp, wk, MPFR_RNDN); 2572 mpfr_abs(tmp, tmp, MPFR_RNDN); 2573 } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor * digits); /* Only need to evaluate sum until weights are < 32 digits of precision */ 2574 mpfr_sub(tmp, sum, osum, MPFR_RNDN); 2575 mpfr_abs(tmp, tmp, MPFR_RNDN); 2576 mpfr_log10(tmp, tmp, MPFR_RNDN); 2577 d1 = mpfr_get_d(tmp, MPFR_RNDN); 2578 mpfr_sub(tmp, sum, psum, MPFR_RNDN); 2579 mpfr_abs(tmp, tmp, MPFR_RNDN); 2580 mpfr_log10(tmp, tmp, MPFR_RNDN); 2581 d2 = mpfr_get_d(tmp, MPFR_RNDN); 2582 mpfr_log10(tmp, maxTerm, MPFR_RNDN); 2583 d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits; 2584 mpfr_log10(tmp, curTerm, MPFR_RNDN); 2585 d4 = mpfr_get_d(tmp, MPFR_RNDN); 2586 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4))); 2587 } while (d < digits && l < 8); 2588 *sol = mpfr_get_d(sum, MPFR_RNDN); 2589 /* Cleanup */ 2590 mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL); 2591 PetscFunctionReturn(PETSC_SUCCESS); 2592 } 2593 #else 2594 2595 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2596 { 2597 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp"); 2598 } 2599 #endif 2600 2601 /*@ 2602 PetscDTTensorQuadratureCreate - create the tensor product quadrature from two lower-dimensional quadratures 2603 2604 Not Collective 2605 2606 Input Parameters: 2607 + q1 - The first quadrature 2608 - q2 - The second quadrature 2609 2610 Output Parameter: 2611 . q - A `PetscQuadrature` object 2612 2613 Level: intermediate 2614 2615 .seealso: `PetscQuadrature`, `PetscDTGaussTensorQuadrature()` 2616 @*/ 2617 PetscErrorCode PetscDTTensorQuadratureCreate(PetscQuadrature q1, PetscQuadrature q2, PetscQuadrature *q) 2618 { 2619 DMPolytopeType ct1, ct2, ct; 2620 const PetscReal *x1, *w1, *x2, *w2; 2621 PetscReal *x, *w; 2622 PetscInt dim1, Nc1, Np1, order1, qa, d1; 2623 PetscInt dim2, Nc2, Np2, order2, qb, d2; 2624 PetscInt dim, Nc, Np, order, qc, d; 2625 2626 PetscFunctionBegin; 2627 PetscValidHeaderSpecific(q1, PETSCQUADRATURE_CLASSID, 1); 2628 PetscValidHeaderSpecific(q2, PETSCQUADRATURE_CLASSID, 2); 2629 PetscAssertPointer(q, 3); 2630 2631 PetscCall(PetscQuadratureGetOrder(q1, &order1)); 2632 PetscCall(PetscQuadratureGetOrder(q2, &order2)); 2633 PetscCheck(order1 == order2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Order1 %" PetscInt_FMT " != %" PetscInt_FMT " Order2", order1, order2); 2634 PetscCall(PetscQuadratureGetData(q1, &dim1, &Nc1, &Np1, &x1, &w1)); 2635 PetscCall(PetscQuadratureGetCellType(q1, &ct1)); 2636 PetscCall(PetscQuadratureGetData(q2, &dim2, &Nc2, &Np2, &x2, &w2)); 2637 PetscCall(PetscQuadratureGetCellType(q2, &ct2)); 2638 PetscCheck(Nc1 == Nc2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "NumComp1 %" PetscInt_FMT " != %" PetscInt_FMT " NumComp2", Nc1, Nc2); 2639 2640 switch (ct1) { 2641 case DM_POLYTOPE_POINT: 2642 ct = ct2; 2643 break; 2644 case DM_POLYTOPE_SEGMENT: 2645 switch (ct2) { 2646 case DM_POLYTOPE_POINT: 2647 ct = ct1; 2648 break; 2649 case DM_POLYTOPE_SEGMENT: 2650 ct = DM_POLYTOPE_QUADRILATERAL; 2651 break; 2652 case DM_POLYTOPE_TRIANGLE: 2653 ct = DM_POLYTOPE_TRI_PRISM; 2654 break; 2655 case DM_POLYTOPE_QUADRILATERAL: 2656 ct = DM_POLYTOPE_HEXAHEDRON; 2657 break; 2658 case DM_POLYTOPE_TETRAHEDRON: 2659 ct = DM_POLYTOPE_UNKNOWN; 2660 break; 2661 case DM_POLYTOPE_HEXAHEDRON: 2662 ct = DM_POLYTOPE_UNKNOWN; 2663 break; 2664 default: 2665 ct = DM_POLYTOPE_UNKNOWN; 2666 } 2667 break; 2668 case DM_POLYTOPE_TRIANGLE: 2669 switch (ct2) { 2670 case DM_POLYTOPE_POINT: 2671 ct = ct1; 2672 break; 2673 case DM_POLYTOPE_SEGMENT: 2674 ct = DM_POLYTOPE_TRI_PRISM; 2675 break; 2676 case DM_POLYTOPE_TRIANGLE: 2677 ct = DM_POLYTOPE_UNKNOWN; 2678 break; 2679 case DM_POLYTOPE_QUADRILATERAL: 2680 ct = DM_POLYTOPE_UNKNOWN; 2681 break; 2682 case DM_POLYTOPE_TETRAHEDRON: 2683 ct = DM_POLYTOPE_UNKNOWN; 2684 break; 2685 case DM_POLYTOPE_HEXAHEDRON: 2686 ct = DM_POLYTOPE_UNKNOWN; 2687 break; 2688 default: 2689 ct = DM_POLYTOPE_UNKNOWN; 2690 } 2691 break; 2692 case DM_POLYTOPE_QUADRILATERAL: 2693 switch (ct2) { 2694 case DM_POLYTOPE_POINT: 2695 ct = ct1; 2696 break; 2697 case DM_POLYTOPE_SEGMENT: 2698 ct = DM_POLYTOPE_HEXAHEDRON; 2699 break; 2700 case DM_POLYTOPE_TRIANGLE: 2701 ct = DM_POLYTOPE_UNKNOWN; 2702 break; 2703 case DM_POLYTOPE_QUADRILATERAL: 2704 ct = DM_POLYTOPE_UNKNOWN; 2705 break; 2706 case DM_POLYTOPE_TETRAHEDRON: 2707 ct = DM_POLYTOPE_UNKNOWN; 2708 break; 2709 case DM_POLYTOPE_HEXAHEDRON: 2710 ct = DM_POLYTOPE_UNKNOWN; 2711 break; 2712 default: 2713 ct = DM_POLYTOPE_UNKNOWN; 2714 } 2715 break; 2716 case DM_POLYTOPE_TETRAHEDRON: 2717 switch (ct2) { 2718 case DM_POLYTOPE_POINT: 2719 ct = ct1; 2720 break; 2721 case DM_POLYTOPE_SEGMENT: 2722 ct = DM_POLYTOPE_UNKNOWN; 2723 break; 2724 case DM_POLYTOPE_TRIANGLE: 2725 ct = DM_POLYTOPE_UNKNOWN; 2726 break; 2727 case DM_POLYTOPE_QUADRILATERAL: 2728 ct = DM_POLYTOPE_UNKNOWN; 2729 break; 2730 case DM_POLYTOPE_TETRAHEDRON: 2731 ct = DM_POLYTOPE_UNKNOWN; 2732 break; 2733 case DM_POLYTOPE_HEXAHEDRON: 2734 ct = DM_POLYTOPE_UNKNOWN; 2735 break; 2736 default: 2737 ct = DM_POLYTOPE_UNKNOWN; 2738 } 2739 break; 2740 case DM_POLYTOPE_HEXAHEDRON: 2741 switch (ct2) { 2742 case DM_POLYTOPE_POINT: 2743 ct = ct1; 2744 break; 2745 case DM_POLYTOPE_SEGMENT: 2746 ct = DM_POLYTOPE_UNKNOWN; 2747 break; 2748 case DM_POLYTOPE_TRIANGLE: 2749 ct = DM_POLYTOPE_UNKNOWN; 2750 break; 2751 case DM_POLYTOPE_QUADRILATERAL: 2752 ct = DM_POLYTOPE_UNKNOWN; 2753 break; 2754 case DM_POLYTOPE_TETRAHEDRON: 2755 ct = DM_POLYTOPE_UNKNOWN; 2756 break; 2757 case DM_POLYTOPE_HEXAHEDRON: 2758 ct = DM_POLYTOPE_UNKNOWN; 2759 break; 2760 default: 2761 ct = DM_POLYTOPE_UNKNOWN; 2762 } 2763 break; 2764 default: 2765 ct = DM_POLYTOPE_UNKNOWN; 2766 } 2767 dim = dim1 + dim2; 2768 Nc = Nc1; 2769 Np = Np1 * Np2; 2770 order = order1; 2771 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 2772 PetscCall(PetscQuadratureSetCellType(*q, ct)); 2773 PetscCall(PetscQuadratureSetOrder(*q, order)); 2774 PetscCall(PetscMalloc1(Np * dim, &x)); 2775 PetscCall(PetscMalloc1(Np, &w)); 2776 for (qa = 0, qc = 0; qa < Np1; ++qa) { 2777 for (qb = 0; qb < Np2; ++qb, ++qc) { 2778 for (d1 = 0, d = 0; d1 < dim1; ++d1, ++d) x[qc * dim + d] = x1[qa * dim1 + d1]; 2779 for (d2 = 0; d2 < dim2; ++d2, ++d) x[qc * dim + d] = x2[qb * dim2 + d2]; 2780 w[qc] = w1[qa] * w2[qb]; 2781 } 2782 } 2783 PetscCall(PetscQuadratureSetData(*q, dim, Nc, Np, x, w)); 2784 PetscFunctionReturn(PETSC_SUCCESS); 2785 } 2786 2787 /* Overwrites A. Can only handle full-rank problems with m>=n 2788 A in column-major format 2789 Ainv in row-major format 2790 tau has length m 2791 worksize must be >= max(1,n) 2792 */ 2793 static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m, PetscInt mstride, PetscInt n, PetscReal *A_in, PetscReal *Ainv_out, PetscScalar *tau, PetscInt worksize, PetscScalar *work) 2794 { 2795 PetscBLASInt M, N, K, lda, ldb, ldwork, info; 2796 PetscScalar *A, *Ainv, *R, *Q, Alpha; 2797 2798 PetscFunctionBegin; 2799 #if defined(PETSC_USE_COMPLEX) 2800 { 2801 PetscInt i, j; 2802 PetscCall(PetscMalloc2(m * n, &A, m * n, &Ainv)); 2803 for (j = 0; j < n; j++) { 2804 for (i = 0; i < m; i++) A[i + m * j] = A_in[i + mstride * j]; 2805 } 2806 mstride = m; 2807 } 2808 #else 2809 A = A_in; 2810 Ainv = Ainv_out; 2811 #endif 2812 2813 PetscCall(PetscBLASIntCast(m, &M)); 2814 PetscCall(PetscBLASIntCast(n, &N)); 2815 PetscCall(PetscBLASIntCast(mstride, &lda)); 2816 PetscCall(PetscBLASIntCast(worksize, &ldwork)); 2817 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 2818 PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info)); 2819 PetscCall(PetscFPTrapPop()); 2820 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error"); 2821 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */ 2822 2823 /* Extract an explicit representation of Q */ 2824 Q = Ainv; 2825 PetscCall(PetscArraycpy(Q, A, mstride * n)); 2826 K = N; /* full rank */ 2827 PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info)); 2828 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error"); 2829 2830 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */ 2831 Alpha = 1.0; 2832 ldb = lda; 2833 PetscCallBLAS("BLAStrsm", BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb)); 2834 /* Ainv is Q, overwritten with inverse */ 2835 2836 #if defined(PETSC_USE_COMPLEX) 2837 { 2838 PetscInt i; 2839 for (i = 0; i < m * n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]); 2840 PetscCall(PetscFree2(A, Ainv)); 2841 } 2842 #endif 2843 PetscFunctionReturn(PETSC_SUCCESS); 2844 } 2845 2846 /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */ 2847 static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval, const PetscReal *x, PetscInt ndegree, const PetscInt *degrees, PetscBool Transpose, PetscReal *B) 2848 { 2849 PetscReal *Bv; 2850 PetscInt i, j; 2851 2852 PetscFunctionBegin; 2853 PetscCall(PetscMalloc1((ninterval + 1) * ndegree, &Bv)); 2854 /* Point evaluation of L_p on all the source vertices */ 2855 PetscCall(PetscDTLegendreEval(ninterval + 1, x, ndegree, degrees, Bv, NULL, NULL)); 2856 /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */ 2857 for (i = 0; i < ninterval; i++) { 2858 for (j = 0; j < ndegree; j++) { 2859 if (Transpose) B[i + ninterval * j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j]; 2860 else B[i * ndegree + j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j]; 2861 } 2862 } 2863 PetscCall(PetscFree(Bv)); 2864 PetscFunctionReturn(PETSC_SUCCESS); 2865 } 2866 2867 /*@ 2868 PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals 2869 2870 Not Collective 2871 2872 Input Parameters: 2873 + degree - degree of reconstruction polynomial 2874 . nsource - number of source intervals 2875 . sourcex - sorted coordinates of source cell boundaries (length `nsource`+1) 2876 . ntarget - number of target intervals 2877 - targetx - sorted coordinates of target cell boundaries (length `ntarget`+1) 2878 2879 Output Parameter: 2880 . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s] 2881 2882 Level: advanced 2883 2884 .seealso: `PetscDTLegendreEval()` 2885 @*/ 2886 PetscErrorCode PetscDTReconstructPoly(PetscInt degree, PetscInt nsource, const PetscReal sourcex[], PetscInt ntarget, const PetscReal targetx[], PetscReal R[]) 2887 { 2888 PetscInt i, j, k, *bdegrees, worksize; 2889 PetscReal xmin, xmax, center, hscale, *sourcey, *targety, *Bsource, *Bsinv, *Btarget; 2890 PetscScalar *tau, *work; 2891 2892 PetscFunctionBegin; 2893 PetscAssertPointer(sourcex, 3); 2894 PetscAssertPointer(targetx, 5); 2895 PetscAssertPointer(R, 6); 2896 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); 2897 if (PetscDefined(USE_DEBUG)) { 2898 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]); 2899 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]); 2900 } 2901 xmin = PetscMin(sourcex[0], targetx[0]); 2902 xmax = PetscMax(sourcex[nsource], targetx[ntarget]); 2903 center = (xmin + xmax) / 2; 2904 hscale = (xmax - xmin) / 2; 2905 worksize = nsource; 2906 PetscCall(PetscMalloc4(degree + 1, &bdegrees, nsource + 1, &sourcey, nsource * (degree + 1), &Bsource, worksize, &work)); 2907 PetscCall(PetscMalloc4(nsource, &tau, nsource * (degree + 1), &Bsinv, ntarget + 1, &targety, ntarget * (degree + 1), &Btarget)); 2908 for (i = 0; i <= nsource; i++) sourcey[i] = (sourcex[i] - center) / hscale; 2909 for (i = 0; i <= degree; i++) bdegrees[i] = i + 1; 2910 PetscCall(PetscDTLegendreIntegrate(nsource, sourcey, degree + 1, bdegrees, PETSC_TRUE, Bsource)); 2911 PetscCall(PetscDTPseudoInverseQR(nsource, nsource, degree + 1, Bsource, Bsinv, tau, nsource, work)); 2912 for (i = 0; i <= ntarget; i++) targety[i] = (targetx[i] - center) / hscale; 2913 PetscCall(PetscDTLegendreIntegrate(ntarget, targety, degree + 1, bdegrees, PETSC_FALSE, Btarget)); 2914 for (i = 0; i < ntarget; i++) { 2915 PetscReal rowsum = 0; 2916 for (j = 0; j < nsource; j++) { 2917 PetscReal sum = 0; 2918 for (k = 0; k < degree + 1; k++) sum += Btarget[i * (degree + 1) + k] * Bsinv[k * nsource + j]; 2919 R[i * nsource + j] = sum; 2920 rowsum += sum; 2921 } 2922 for (j = 0; j < nsource; j++) R[i * nsource + j] /= rowsum; /* normalize each row */ 2923 } 2924 PetscCall(PetscFree4(bdegrees, sourcey, Bsource, work)); 2925 PetscCall(PetscFree4(tau, Bsinv, targety, Btarget)); 2926 PetscFunctionReturn(PETSC_SUCCESS); 2927 } 2928 2929 /*@ 2930 PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points 2931 2932 Not Collective 2933 2934 Input Parameters: 2935 + n - the number of GLL nodes 2936 . nodes - the GLL nodes 2937 . weights - the GLL weights 2938 - f - the function values at the nodes 2939 2940 Output Parameter: 2941 . in - the value of the integral 2942 2943 Level: beginner 2944 2945 .seealso: `PetscDTGaussLobattoLegendreQuadrature()` 2946 @*/ 2947 PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n, PetscReal nodes[], PetscReal weights[], const PetscReal f[], PetscReal *in) 2948 { 2949 PetscInt i; 2950 2951 PetscFunctionBegin; 2952 *in = 0.; 2953 for (i = 0; i < n; i++) *in += f[i] * f[i] * weights[i]; 2954 PetscFunctionReturn(PETSC_SUCCESS); 2955 } 2956 2957 /*@C 2958 PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element 2959 2960 Not Collective 2961 2962 Input Parameters: 2963 + n - the number of GLL nodes 2964 . nodes - the GLL nodes, of length `n` 2965 - weights - the GLL weights, of length `n` 2966 2967 Output Parameter: 2968 . AA - the stiffness element, of size `n` by `n` 2969 2970 Level: beginner 2971 2972 Notes: 2973 Destroy this with `PetscGaussLobattoLegendreElementLaplacianDestroy()` 2974 2975 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) 2976 2977 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()` 2978 @*/ 2979 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 2980 { 2981 PetscReal **A; 2982 const PetscReal *gllnodes = nodes; 2983 const PetscInt p = n - 1; 2984 PetscReal z0, z1, z2 = -1, x, Lpj, Lpr; 2985 PetscInt i, j, nn, r; 2986 2987 PetscFunctionBegin; 2988 PetscCall(PetscMalloc1(n, &A)); 2989 PetscCall(PetscMalloc1(n * n, &A[0])); 2990 for (i = 1; i < n; i++) A[i] = A[i - 1] + n; 2991 2992 for (j = 1; j < p; j++) { 2993 x = gllnodes[j]; 2994 z0 = 1.; 2995 z1 = x; 2996 for (nn = 1; nn < p; nn++) { 2997 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 2998 z0 = z1; 2999 z1 = z2; 3000 } 3001 Lpj = z2; 3002 for (r = 1; r < p; r++) { 3003 if (r == j) { 3004 A[j][j] = 2. / (3. * (1. - gllnodes[j] * gllnodes[j]) * Lpj * Lpj); 3005 } else { 3006 x = gllnodes[r]; 3007 z0 = 1.; 3008 z1 = x; 3009 for (nn = 1; nn < p; nn++) { 3010 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 3011 z0 = z1; 3012 z1 = z2; 3013 } 3014 Lpr = z2; 3015 A[r][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * Lpr * (gllnodes[j] - gllnodes[r]) * (gllnodes[j] - gllnodes[r])); 3016 } 3017 } 3018 } 3019 for (j = 1; j < p + 1; j++) { 3020 x = gllnodes[j]; 3021 z0 = 1.; 3022 z1 = x; 3023 for (nn = 1; nn < p; nn++) { 3024 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 3025 z0 = z1; 3026 z1 = z2; 3027 } 3028 Lpj = z2; 3029 A[j][0] = 4. * PetscPowRealInt(-1., p) / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. + gllnodes[j]) * (1. + gllnodes[j])); 3030 A[0][j] = A[j][0]; 3031 } 3032 for (j = 0; j < p; j++) { 3033 x = gllnodes[j]; 3034 z0 = 1.; 3035 z1 = x; 3036 for (nn = 1; nn < p; nn++) { 3037 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 3038 z0 = z1; 3039 z1 = z2; 3040 } 3041 Lpj = z2; 3042 3043 A[p][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. - gllnodes[j]) * (1. - gllnodes[j])); 3044 A[j][p] = A[p][j]; 3045 } 3046 A[0][0] = 0.5 + (((PetscReal)p) * (((PetscReal)p) + 1.) - 2.) / 6.; 3047 A[p][p] = A[0][0]; 3048 *AA = A; 3049 PetscFunctionReturn(PETSC_SUCCESS); 3050 } 3051 3052 /*@C 3053 PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element created with `PetscGaussLobattoLegendreElementLaplacianCreate()` 3054 3055 Not Collective 3056 3057 Input Parameters: 3058 + n - the number of GLL nodes 3059 . nodes - the GLL nodes, ignored 3060 . weights - the GLL weightss, ignored 3061 - AA - the stiffness element from `PetscGaussLobattoLegendreElementLaplacianCreate()` 3062 3063 Level: beginner 3064 3065 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()` 3066 @*/ 3067 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 3068 { 3069 PetscFunctionBegin; 3070 PetscCall(PetscFree((*AA)[0])); 3071 PetscCall(PetscFree(*AA)); 3072 *AA = NULL; 3073 PetscFunctionReturn(PETSC_SUCCESS); 3074 } 3075 3076 /*@C 3077 PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element 3078 3079 Not Collective 3080 3081 Input Parameters: 3082 + n - the number of GLL nodes 3083 . nodes - the GLL nodes, of length `n` 3084 - weights - the GLL weights, of length `n` 3085 3086 Output Parameters: 3087 + AA - the stiffness element, of dimension `n` by `n` 3088 - AAT - the transpose of AA (pass in `NULL` if you do not need this array), of dimension `n` by `n` 3089 3090 Level: beginner 3091 3092 Notes: 3093 Destroy this with `PetscGaussLobattoLegendreElementGradientDestroy()` 3094 3095 You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row-oriented 3096 3097 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`, `PetscGaussLobattoLegendreElementGradientDestroy()` 3098 @*/ 3099 PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA, PetscReal ***AAT) 3100 { 3101 PetscReal **A, **AT = NULL; 3102 const PetscReal *gllnodes = nodes; 3103 const PetscInt p = n - 1; 3104 PetscReal Li, Lj, d0; 3105 PetscInt i, j; 3106 3107 PetscFunctionBegin; 3108 PetscCall(PetscMalloc1(n, &A)); 3109 PetscCall(PetscMalloc1(n * n, &A[0])); 3110 for (i = 1; i < n; i++) A[i] = A[i - 1] + n; 3111 3112 if (AAT) { 3113 PetscCall(PetscMalloc1(n, &AT)); 3114 PetscCall(PetscMalloc1(n * n, &AT[0])); 3115 for (i = 1; i < n; i++) AT[i] = AT[i - 1] + n; 3116 } 3117 3118 if (n == 1) A[0][0] = 0.; 3119 d0 = (PetscReal)p * ((PetscReal)p + 1.) / 4.; 3120 for (i = 0; i < n; i++) { 3121 for (j = 0; j < n; j++) { 3122 A[i][j] = 0.; 3123 PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li)); 3124 PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj)); 3125 if (i != j) A[i][j] = Li / (Lj * (gllnodes[i] - gllnodes[j])); 3126 if ((j == i) && (i == 0)) A[i][j] = -d0; 3127 if (j == i && i == p) A[i][j] = d0; 3128 if (AT) AT[j][i] = A[i][j]; 3129 } 3130 } 3131 if (AAT) *AAT = AT; 3132 *AA = A; 3133 PetscFunctionReturn(PETSC_SUCCESS); 3134 } 3135 3136 /*@C 3137 PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with `PetscGaussLobattoLegendreElementGradientCreate()` 3138 3139 Not Collective 3140 3141 Input Parameters: 3142 + n - the number of GLL nodes 3143 . nodes - the GLL nodes, ignored 3144 . weights - the GLL weights, ignored 3145 . AA - the stiffness element obtained with `PetscGaussLobattoLegendreElementGradientCreate()` 3146 - AAT - the transpose of the element obtained with `PetscGaussLobattoLegendreElementGradientCreate()` 3147 3148 Level: beginner 3149 3150 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionCreate()` 3151 @*/ 3152 PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA, PetscReal ***AAT) 3153 { 3154 PetscFunctionBegin; 3155 PetscCall(PetscFree((*AA)[0])); 3156 PetscCall(PetscFree(*AA)); 3157 *AA = NULL; 3158 if (AAT) { 3159 PetscCall(PetscFree((*AAT)[0])); 3160 PetscCall(PetscFree(*AAT)); 3161 *AAT = NULL; 3162 } 3163 PetscFunctionReturn(PETSC_SUCCESS); 3164 } 3165 3166 /*@C 3167 PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element 3168 3169 Not Collective 3170 3171 Input Parameters: 3172 + n - the number of GLL nodes 3173 . nodes - the GLL nodes, of length `n` 3174 - weights - the GLL weights, of length `n` 3175 3176 Output Parameter: 3177 . AA - the stiffness element, of dimension `n` by `n` 3178 3179 Level: beginner 3180 3181 Notes: 3182 Destroy this with `PetscGaussLobattoLegendreElementAdvectionDestroy()` 3183 3184 This is the same as the Gradient operator multiplied by the diagonal mass matrix 3185 3186 You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row-oriented 3187 3188 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionDestroy()` 3189 @*/ 3190 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 3191 { 3192 PetscReal **D; 3193 const PetscReal *gllweights = weights; 3194 const PetscInt glln = n; 3195 PetscInt i, j; 3196 3197 PetscFunctionBegin; 3198 PetscCall(PetscGaussLobattoLegendreElementGradientCreate(n, nodes, weights, &D, NULL)); 3199 for (i = 0; i < glln; i++) { 3200 for (j = 0; j < glln; j++) D[i][j] = gllweights[i] * D[i][j]; 3201 } 3202 *AA = D; 3203 PetscFunctionReturn(PETSC_SUCCESS); 3204 } 3205 3206 /*@C 3207 PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element created with `PetscGaussLobattoLegendreElementAdvectionCreate()` 3208 3209 Not Collective 3210 3211 Input Parameters: 3212 + n - the number of GLL nodes 3213 . nodes - the GLL nodes, ignored 3214 . weights - the GLL weights, ignored 3215 - AA - advection obtained with `PetscGaussLobattoLegendreElementAdvectionCreate()` 3216 3217 Level: beginner 3218 3219 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementAdvectionCreate()` 3220 @*/ 3221 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 3222 { 3223 PetscFunctionBegin; 3224 PetscCall(PetscFree((*AA)[0])); 3225 PetscCall(PetscFree(*AA)); 3226 *AA = NULL; 3227 PetscFunctionReturn(PETSC_SUCCESS); 3228 } 3229 3230 PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 3231 { 3232 PetscReal **A; 3233 const PetscReal *gllweights = weights; 3234 const PetscInt glln = n; 3235 PetscInt i, j; 3236 3237 PetscFunctionBegin; 3238 PetscCall(PetscMalloc1(glln, &A)); 3239 PetscCall(PetscMalloc1(glln * glln, &A[0])); 3240 for (i = 1; i < glln; i++) A[i] = A[i - 1] + glln; 3241 if (glln == 1) A[0][0] = 0.; 3242 for (i = 0; i < glln; i++) { 3243 for (j = 0; j < glln; j++) { 3244 A[i][j] = 0.; 3245 if (j == i) A[i][j] = gllweights[i]; 3246 } 3247 } 3248 *AA = A; 3249 PetscFunctionReturn(PETSC_SUCCESS); 3250 } 3251 3252 PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 3253 { 3254 PetscFunctionBegin; 3255 PetscCall(PetscFree((*AA)[0])); 3256 PetscCall(PetscFree(*AA)); 3257 *AA = NULL; 3258 PetscFunctionReturn(PETSC_SUCCESS); 3259 } 3260 3261 /*@ 3262 PetscDTIndexToBary - convert an index into a barycentric coordinate. 3263 3264 Input Parameters: 3265 + 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) 3266 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 3267 - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum) 3268 3269 Output Parameter: 3270 . coord - will be filled with the barycentric coordinate, of length `n` 3271 3272 Level: beginner 3273 3274 Note: 3275 The indices map to barycentric coordinates in lexicographic order, where the first index is the 3276 least significant and the last index is the most significant. 3277 3278 .seealso: `PetscDTBaryToIndex()` 3279 @*/ 3280 PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[]) 3281 { 3282 PetscInt c, d, s, total, subtotal, nexttotal; 3283 3284 PetscFunctionBeginHot; 3285 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 3286 PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative"); 3287 if (!len) { 3288 if (!sum && !index) PetscFunctionReturn(PETSC_SUCCESS); 3289 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 3290 } 3291 for (c = 1, total = 1; c <= len; c++) { 3292 /* total is the number of ways to have a tuple of length c with sum */ 3293 if (index < total) break; 3294 total = (total * (sum + c)) / c; 3295 } 3296 PetscCheck(c <= len, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range"); 3297 for (d = c; d < len; d++) coord[d] = 0; 3298 for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) { 3299 /* subtotal is the number of ways to have a tuple of length c with sum s */ 3300 /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */ 3301 if ((index + subtotal) >= total) { 3302 coord[--c] = sum - s; 3303 index -= (total - subtotal); 3304 sum = s; 3305 total = nexttotal; 3306 subtotal = 1; 3307 nexttotal = 1; 3308 s = 0; 3309 } else { 3310 subtotal = (subtotal * (c + s)) / (s + 1); 3311 nexttotal = (nexttotal * (c - 1 + s)) / (s + 1); 3312 s++; 3313 } 3314 } 3315 PetscFunctionReturn(PETSC_SUCCESS); 3316 } 3317 3318 /*@ 3319 PetscDTBaryToIndex - convert a barycentric coordinate to an index 3320 3321 Input Parameters: 3322 + 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) 3323 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 3324 - coord - a barycentric coordinate with the given length `len` and `sum` 3325 3326 Output Parameter: 3327 . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum) 3328 3329 Level: beginner 3330 3331 Note: 3332 The indices map to barycentric coordinates in lexicographic order, where the first index is the 3333 least significant and the last index is the most significant. 3334 3335 .seealso: `PetscDTIndexToBary` 3336 @*/ 3337 PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index) 3338 { 3339 PetscInt c; 3340 PetscInt i; 3341 PetscInt total; 3342 3343 PetscFunctionBeginHot; 3344 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 3345 if (!len) { 3346 if (!sum) { 3347 *index = 0; 3348 PetscFunctionReturn(PETSC_SUCCESS); 3349 } 3350 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 3351 } 3352 for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c; 3353 i = total - 1; 3354 c = len - 1; 3355 sum -= coord[c]; 3356 while (sum > 0) { 3357 PetscInt subtotal; 3358 PetscInt s; 3359 3360 for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s; 3361 i -= subtotal; 3362 sum -= coord[--c]; 3363 } 3364 *index = i; 3365 PetscFunctionReturn(PETSC_SUCCESS); 3366 } 3367 3368 /*@ 3369 PetscQuadratureComputePermutations - Compute permutations of quadrature points corresponding to domain orientations 3370 3371 Input Parameter: 3372 . quad - The `PetscQuadrature` 3373 3374 Output Parameters: 3375 + Np - The number of domain orientations 3376 - perm - An array of `IS` permutations, one for ech orientation, 3377 3378 Level: developer 3379 3380 .seealso: `PetscQuadratureSetCellType()`, `PetscQuadrature` 3381 @*/ 3382 PetscErrorCode PetscQuadratureComputePermutations(PetscQuadrature quad, PeOp PetscInt *Np, IS *perm[]) 3383 { 3384 DMPolytopeType ct; 3385 const PetscReal *xq, *wq; 3386 PetscInt dim, qdim, d, Na, o, Nq, q, qp; 3387 3388 PetscFunctionBegin; 3389 PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &xq, &wq)); 3390 PetscCall(PetscQuadratureGetCellType(quad, &ct)); 3391 dim = DMPolytopeTypeGetDim(ct); 3392 Na = DMPolytopeTypeGetNumArrangements(ct); 3393 PetscCall(PetscMalloc1(Na, perm)); 3394 if (Np) *Np = Na; 3395 Na /= 2; 3396 for (o = -Na; o < Na; ++o) { 3397 DM refdm; 3398 PetscInt *idx; 3399 PetscReal xi0[3] = {-1., -1., -1.}, v0[3], J[9], detJ, txq[3]; 3400 PetscBool flg; 3401 3402 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm)); 3403 PetscCall(DMPlexOrientPoint(refdm, 0, o)); 3404 PetscCall(DMPlexComputeCellGeometryFEM(refdm, 0, NULL, v0, J, NULL, &detJ)); 3405 PetscCall(PetscMalloc1(Nq, &idx)); 3406 for (q = 0; q < Nq; ++q) { 3407 CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q * dim], txq); 3408 for (qp = 0; qp < Nq; ++qp) { 3409 PetscReal diff = 0.; 3410 3411 for (d = 0; d < dim; ++d) diff += PetscAbsReal(txq[d] - xq[qp * dim + d]); 3412 if (diff < PETSC_SMALL) break; 3413 } 3414 PetscCheck(qp < Nq, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Transformed quad point %" PetscInt_FMT " does not match another quad point", q); 3415 idx[q] = qp; 3416 } 3417 PetscCall(DMDestroy(&refdm)); 3418 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nq, idx, PETSC_OWN_POINTER, &(*perm)[o + Na])); 3419 PetscCall(ISGetInfo((*perm)[o + Na], IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg)); 3420 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ordering for orientation %" PetscInt_FMT " was not a permutation", o); 3421 PetscCall(ISSetPermutation((*perm)[o + Na])); 3422 } 3423 if (!Na) (*perm)[0] = NULL; 3424 PetscFunctionReturn(PETSC_SUCCESS); 3425 } 3426 3427 /*@ 3428 PetscDTCreateQuadratureByCell - Create default quadrature for a given cell 3429 3430 Not collective 3431 3432 Input Parameters: 3433 + ct - The integration domain 3434 . qorder - The desired quadrature order 3435 - qtype - The type of simplex quadrature, or PETSCDTSIMPLEXQUAD_DEFAULT 3436 3437 Output Parameters: 3438 + q - The cell quadrature 3439 - fq - The face quadrature 3440 3441 Level: developer 3442 3443 .seealso: `PetscDTCreateDefaultQuadrature()`, `PetscFECreateDefault()`, `PetscDTGaussTensorQuadrature()`, `PetscDTSimplexQuadrature()`, `PetscDTTensorQuadratureCreate()` 3444 @*/ 3445 PetscErrorCode PetscDTCreateQuadratureByCell(DMPolytopeType ct, PetscInt qorder, PetscDTSimplexQuadratureType qtype, PetscQuadrature *q, PetscQuadrature *fq) 3446 { 3447 const PetscInt quadPointsPerEdge = PetscMax(qorder + 1, 1); 3448 const PetscInt dim = DMPolytopeTypeGetDim(ct); 3449 3450 PetscFunctionBegin; 3451 switch (ct) { 3452 case DM_POLYTOPE_SEGMENT: 3453 case DM_POLYTOPE_POINT_PRISM_TENSOR: 3454 case DM_POLYTOPE_QUADRILATERAL: 3455 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3456 case DM_POLYTOPE_HEXAHEDRON: 3457 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3458 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q)); 3459 PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq)); 3460 break; 3461 case DM_POLYTOPE_TRIANGLE: 3462 case DM_POLYTOPE_TETRAHEDRON: 3463 PetscCall(PetscDTSimplexQuadrature(dim, 2 * qorder, qtype, q)); 3464 PetscCall(PetscDTSimplexQuadrature(dim - 1, 2 * qorder, qtype, fq)); 3465 break; 3466 case DM_POLYTOPE_TRI_PRISM: 3467 case DM_POLYTOPE_TRI_PRISM_TENSOR: { 3468 PetscQuadrature q1, q2; 3469 3470 // TODO: this should be able to use symmetric rules, but doing so causes tests to fail 3471 PetscCall(PetscDTSimplexQuadrature(2, 2 * qorder, PETSCDTSIMPLEXQUAD_CONIC, &q1)); 3472 PetscCall(PetscDTGaussTensorQuadrature(1, 1, quadPointsPerEdge, -1.0, 1.0, &q2)); 3473 PetscCall(PetscDTTensorQuadratureCreate(q1, q2, q)); 3474 PetscCall(PetscQuadratureDestroy(&q2)); 3475 *fq = q1; 3476 /* TODO Need separate quadratures for each face */ 3477 } break; 3478 default: 3479 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No quadrature for celltype %s", DMPolytopeTypes[PetscMin(ct, DM_POLYTOPE_UNKNOWN)]); 3480 } 3481 PetscFunctionReturn(PETSC_SUCCESS); 3482 } 3483 3484 /*@ 3485 PetscDTCreateDefaultQuadrature - Create default quadrature for a given cell 3486 3487 Not collective 3488 3489 Input Parameters: 3490 + ct - The integration domain 3491 - qorder - The desired quadrature order 3492 3493 Output Parameters: 3494 + q - The cell quadrature 3495 - fq - The face quadrature 3496 3497 Level: developer 3498 3499 .seealso: `PetscDTCreateQuadratureByCell()`, `PetscFECreateDefault()`, `PetscDTGaussTensorQuadrature()`, `PetscDTSimplexQuadrature()`, `PetscDTTensorQuadratureCreate()` 3500 @*/ 3501 PetscErrorCode PetscDTCreateDefaultQuadrature(DMPolytopeType ct, PetscInt qorder, PetscQuadrature *q, PetscQuadrature *fq) 3502 { 3503 PetscFunctionBegin; 3504 PetscCall(PetscDTCreateQuadratureByCell(ct, qorder, PETSCDTSIMPLEXQUAD_DEFAULT, q, fq)); 3505 PetscFunctionReturn(PETSC_SUCCESS); 3506 } 3507