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 isascii; 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, &isascii)); 637 PetscCall(PetscViewerASCIIPushTab(viewer)); 638 if (isascii) 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 PetscReal binom1, binom2; 1489 PetscInt alphai = (PetscInt)alpha; 1490 PetscInt betai = (PetscInt)beta; 1491 1492 PetscCheck((PetscReal)alphai == alpha && (PetscReal)betai == beta, PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable."); 1493 PetscCall(PetscDTBinomial(m + b, b, &binom1)); 1494 PetscCall(PetscDTBinomial(m + a + b, b, &binom2)); 1495 grb = 1. / (binom1 * binom2); 1496 PetscCall(PetscDTBinomial(m + a, a, &binom1)); 1497 PetscCall(PetscDTBinomial(m + a + b, a, &binom2)); 1498 gra = 1. / (binom1 * binom2); 1499 } 1500 #endif 1501 *leftw = twoab1 * grb / b; 1502 *rightw = twoab1 * gra / a; 1503 PetscFunctionReturn(PETSC_SUCCESS); 1504 } 1505 1506 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x. 1507 Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */ 1508 static inline PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 1509 { 1510 PetscReal pn1, pn2; 1511 PetscReal cnm1, cnm1x, cnm2; 1512 PetscInt k; 1513 1514 PetscFunctionBegin; 1515 if (!n) { 1516 *P = 1.0; 1517 PetscFunctionReturn(PETSC_SUCCESS); 1518 } 1519 PetscDTJacobiRecurrence_Internal(1, a, b, cnm1, cnm1x, cnm2); 1520 pn2 = 1.; 1521 pn1 = cnm1 + cnm1x * x; 1522 if (n == 1) { 1523 *P = pn1; 1524 PetscFunctionReturn(PETSC_SUCCESS); 1525 } 1526 *P = 0.0; 1527 for (k = 2; k < n + 1; ++k) { 1528 PetscDTJacobiRecurrence_Internal(k, a, b, cnm1, cnm1x, cnm2); 1529 1530 *P = (cnm1 + cnm1x * x) * pn1 - cnm2 * pn2; 1531 pn2 = pn1; 1532 pn1 = *P; 1533 } 1534 PetscFunctionReturn(PETSC_SUCCESS); 1535 } 1536 1537 /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */ 1538 static inline PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P) 1539 { 1540 PetscReal nP; 1541 PetscInt i; 1542 1543 PetscFunctionBegin; 1544 *P = 0.0; 1545 if (k > n) PetscFunctionReturn(PETSC_SUCCESS); 1546 PetscCall(PetscDTComputeJacobi(a + k, b + k, n - k, x, &nP)); 1547 for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5; 1548 *P = nP; 1549 PetscFunctionReturn(PETSC_SUCCESS); 1550 } 1551 1552 static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[]) 1553 { 1554 PetscInt maxIter = 100; 1555 PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON)); 1556 PetscReal a1, a6, gf; 1557 PetscInt k; 1558 1559 PetscFunctionBegin; 1560 a1 = PetscPowReal(2.0, a + b + 1); 1561 #if defined(PETSC_HAVE_LGAMMA) 1562 { 1563 PetscReal a2, a3, a4, a5; 1564 a2 = PetscLGamma(a + npoints + 1); 1565 a3 = PetscLGamma(b + npoints + 1); 1566 a4 = PetscLGamma(a + b + npoints + 1); 1567 a5 = PetscLGamma(npoints + 1); 1568 gf = PetscExpReal(a2 + a3 - (a4 + a5)); 1569 } 1570 #else 1571 { 1572 PetscInt ia, ib; 1573 1574 ia = (PetscInt)a; 1575 ib = (PetscInt)b; 1576 gf = 1.; 1577 if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */ 1578 for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k); 1579 } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */ 1580 for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k); 1581 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable."); 1582 } 1583 #endif 1584 1585 a6 = a1 * gf; 1586 /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses. 1587 Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */ 1588 for (k = 0; k < npoints; ++k) { 1589 PetscReal r = PetscCosReal(PETSC_PI * (1. - (4. * k + 3. + 2. * b) / (4. * npoints + 2. * (a + b + 1.)))), dP; 1590 PetscInt j; 1591 1592 if (k > 0) r = 0.5 * (r + x[k - 1]); 1593 for (j = 0; j < maxIter; ++j) { 1594 PetscReal s = 0.0, delta, f, fp; 1595 PetscInt i; 1596 1597 for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]); 1598 PetscCall(PetscDTComputeJacobi(a, b, npoints, r, &f)); 1599 PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp)); 1600 delta = f / (fp - f * s); 1601 r = r - delta; 1602 if (PetscAbsReal(delta) < eps) break; 1603 } 1604 x[k] = r; 1605 PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP)); 1606 w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP); 1607 } 1608 PetscFunctionReturn(PETSC_SUCCESS); 1609 } 1610 1611 /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi 1612 * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */ 1613 static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s) 1614 { 1615 PetscInt i; 1616 1617 PetscFunctionBegin; 1618 for (i = 0; i < nPoints; i++) { 1619 PetscReal A, B, C; 1620 1621 PetscDTJacobiRecurrence_Internal(i + 1, a, b, A, B, C); 1622 d[i] = -A / B; 1623 if (i) s[i - 1] *= C / B; 1624 if (i < nPoints - 1) s[i] = 1. / B; 1625 } 1626 PetscFunctionReturn(PETSC_SUCCESS); 1627 } 1628 1629 static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w) 1630 { 1631 PetscReal mu0; 1632 PetscReal ga, gb, gab; 1633 PetscInt i; 1634 1635 PetscFunctionBegin; 1636 PetscCall(PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite)); 1637 1638 #if defined(PETSC_HAVE_TGAMMA) 1639 ga = PetscTGamma(a + 1); 1640 gb = PetscTGamma(b + 1); 1641 gab = PetscTGamma(a + b + 2); 1642 #else 1643 { 1644 PetscInt ia = (PetscInt)a, ib = (PetscInt)b; 1645 1646 PetscCheck(ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "tgamma() - math routine is unavailable."); 1647 /* All gamma(x) terms are (x-1)! terms */ 1648 PetscCall(PetscDTFactorial(ia, &ga)); 1649 PetscCall(PetscDTFactorial(ib, &gb)); 1650 PetscCall(PetscDTFactorial(ia + ib + 1, &gab)); 1651 } 1652 #endif 1653 mu0 = PetscPowReal(2., a + b + 1.) * ga * gb / gab; 1654 1655 #if defined(PETSCDTGAUSSIANQUADRATURE_EIG) 1656 { 1657 PetscReal *diag, *subdiag; 1658 PetscScalar *V; 1659 1660 PetscCall(PetscMalloc2(npoints, &diag, npoints, &subdiag)); 1661 PetscCall(PetscMalloc1(npoints * npoints, &V)); 1662 PetscCall(PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag)); 1663 for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]); 1664 PetscCall(PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V)); 1665 for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0; 1666 PetscCall(PetscFree(V)); 1667 PetscCall(PetscFree2(diag, subdiag)); 1668 } 1669 #else 1670 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found"); 1671 #endif 1672 { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the 1673 eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that 1674 the eigenvalues are sorted */ 1675 PetscBool sorted; 1676 1677 PetscCall(PetscSortedReal(npoints, x, &sorted)); 1678 if (!sorted) { 1679 PetscInt *order, i; 1680 PetscReal *tmp; 1681 1682 PetscCall(PetscMalloc2(npoints, &order, npoints, &tmp)); 1683 for (i = 0; i < npoints; i++) order[i] = i; 1684 PetscCall(PetscSortRealWithPermutation(npoints, x, order)); 1685 PetscCall(PetscArraycpy(tmp, x, npoints)); 1686 for (i = 0; i < npoints; i++) x[i] = tmp[order[i]]; 1687 PetscCall(PetscArraycpy(tmp, w, npoints)); 1688 for (i = 0; i < npoints; i++) w[i] = tmp[order[i]]; 1689 PetscCall(PetscFree2(order, tmp)); 1690 } 1691 } 1692 PetscFunctionReturn(PETSC_SUCCESS); 1693 } 1694 1695 static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1696 { 1697 PetscFunctionBegin; 1698 PetscCheck(npoints >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive"); 1699 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1700 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1."); 1701 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1."); 1702 1703 if (newton) PetscCall(PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w)); 1704 else PetscCall(PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w)); 1705 if (alpha == beta) { /* symmetrize */ 1706 PetscInt i; 1707 for (i = 0; i < (npoints + 1) / 2; i++) { 1708 PetscInt j = npoints - 1 - i; 1709 PetscReal xi = x[i]; 1710 PetscReal xj = x[j]; 1711 PetscReal wi = w[i]; 1712 PetscReal wj = w[j]; 1713 1714 x[i] = (xi - xj) / 2.; 1715 x[j] = (xj - xi) / 2.; 1716 w[i] = w[j] = (wi + wj) / 2.; 1717 } 1718 } 1719 PetscFunctionReturn(PETSC_SUCCESS); 1720 } 1721 1722 /*@ 1723 PetscDTGaussJacobiQuadrature - quadrature for the interval $[a, b]$ with the weight function 1724 $(x-a)^\alpha (x-b)^\beta$. 1725 1726 Not Collective 1727 1728 Input Parameters: 1729 + npoints - the number of points in the quadrature rule 1730 . a - the left endpoint of the interval 1731 . b - the right endpoint of the interval 1732 . alpha - the left exponent 1733 - beta - the right exponent 1734 1735 Output Parameters: 1736 + x - array of length `npoints`, the locations of the quadrature points 1737 - w - array of length `npoints`, the weights of the quadrature points 1738 1739 Level: intermediate 1740 1741 Note: 1742 This quadrature rule is exact for polynomials up to degree 2*`npoints` - 1. 1743 1744 .seealso: `PetscDTGaussQuadrature()` 1745 @*/ 1746 PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1747 { 1748 PetscInt i; 1749 1750 PetscFunctionBegin; 1751 PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal)); 1752 if (a != -1. || b != 1.) { /* shift */ 1753 for (i = 0; i < npoints; i++) { 1754 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1755 w[i] *= (b - a) / 2.; 1756 } 1757 } 1758 PetscFunctionReturn(PETSC_SUCCESS); 1759 } 1760 1761 static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1762 { 1763 PetscInt i; 1764 1765 PetscFunctionBegin; 1766 PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive"); 1767 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1768 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1."); 1769 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1."); 1770 1771 x[0] = -1.; 1772 x[npoints - 1] = 1.; 1773 if (npoints > 2) PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints - 2, alpha + 1., beta + 1., &x[1], &w[1], newton)); 1774 for (i = 1; i < npoints - 1; i++) w[i] /= (1. - x[i] * x[i]); 1775 PetscCall(PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints - 1])); 1776 PetscFunctionReturn(PETSC_SUCCESS); 1777 } 1778 1779 /*@ 1780 PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval $[a, b]$ with the weight function 1781 $(x-a)^\alpha (x-b)^\beta$, with endpoints `a` and `b` included as quadrature points. 1782 1783 Not Collective 1784 1785 Input Parameters: 1786 + npoints - the number of points in the quadrature rule 1787 . a - the left endpoint of the interval 1788 . b - the right endpoint of the interval 1789 . alpha - the left exponent 1790 - beta - the right exponent 1791 1792 Output Parameters: 1793 + x - array of length `npoints`, the locations of the quadrature points 1794 - w - array of length `npoints`, the weights of the quadrature points 1795 1796 Level: intermediate 1797 1798 Note: 1799 This quadrature rule is exact for polynomials up to degree 2*`npoints` - 3. 1800 1801 .seealso: `PetscDTGaussJacobiQuadrature()` 1802 @*/ 1803 PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1804 { 1805 PetscInt i; 1806 1807 PetscFunctionBegin; 1808 PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal)); 1809 if (a != -1. || b != 1.) { /* shift */ 1810 for (i = 0; i < npoints; i++) { 1811 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1812 w[i] *= (b - a) / 2.; 1813 } 1814 } 1815 PetscFunctionReturn(PETSC_SUCCESS); 1816 } 1817 1818 /*@ 1819 PetscDTGaussQuadrature - create Gauss-Legendre quadrature 1820 1821 Not Collective 1822 1823 Input Parameters: 1824 + npoints - number of points 1825 . a - left end of interval (often-1) 1826 - b - right end of interval (often +1) 1827 1828 Output Parameters: 1829 + x - quadrature points 1830 - w - quadrature weights 1831 1832 Level: intermediate 1833 1834 Note: 1835 See {cite}`golub1969calculation` 1836 1837 .seealso: `PetscDTLegendreEval()`, `PetscDTGaussJacobiQuadrature()` 1838 @*/ 1839 PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[]) 1840 { 1841 PetscInt i; 1842 1843 PetscFunctionBegin; 1844 PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal)); 1845 if (a != -1. || b != 1.) { /* shift */ 1846 for (i = 0; i < npoints; i++) { 1847 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1848 w[i] *= (b - a) / 2.; 1849 } 1850 } 1851 PetscFunctionReturn(PETSC_SUCCESS); 1852 } 1853 1854 /*@ 1855 PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre 1856 nodes of a given size on the domain $[-1,1]$ 1857 1858 Not Collective 1859 1860 Input Parameters: 1861 + npoints - number of grid nodes 1862 - type - `PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA` or `PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON` 1863 1864 Output Parameters: 1865 + x - quadrature points, pass in an array of length `npoints` 1866 - w - quadrature weights, pass in an array of length `npoints` 1867 1868 Level: intermediate 1869 1870 Notes: 1871 For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not 1872 close enough to the desired solution 1873 1874 These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes 1875 1876 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 1877 1878 .seealso: `PetscDTGaussQuadrature()`, `PetscGaussLobattoLegendreCreateType` 1879 1880 @*/ 1881 PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints, PetscGaussLobattoLegendreCreateType type, PetscReal x[], PetscReal w[]) 1882 { 1883 PetscBool newton; 1884 1885 PetscFunctionBegin; 1886 PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide at least 2 grid points per element"); 1887 newton = (PetscBool)(type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON); 1888 PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton)); 1889 PetscFunctionReturn(PETSC_SUCCESS); 1890 } 1891 1892 /*@ 1893 PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature 1894 1895 Not Collective 1896 1897 Input Parameters: 1898 + dim - The spatial dimension 1899 . Nc - The number of components 1900 . npoints - number of points in one dimension 1901 . a - left end of interval (often-1) 1902 - b - right end of interval (often +1) 1903 1904 Output Parameter: 1905 . q - A `PetscQuadrature` object 1906 1907 Level: intermediate 1908 1909 .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()` 1910 @*/ 1911 PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 1912 { 1913 DMPolytopeType ct; 1914 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints * PetscSqr(npoints) : PetscSqr(npoints) : npoints; 1915 PetscReal *x, *w, *xw, *ww; 1916 1917 PetscFunctionBegin; 1918 PetscCall(PetscMalloc1(totpoints * dim, &x)); 1919 PetscCall(PetscMalloc1(totpoints * Nc, &w)); 1920 /* Set up the Golub-Welsch system */ 1921 switch (dim) { 1922 case 0: 1923 ct = DM_POLYTOPE_POINT; 1924 PetscCall(PetscFree(x)); 1925 PetscCall(PetscFree(w)); 1926 PetscCall(PetscMalloc1(1, &x)); 1927 PetscCall(PetscMalloc1(Nc, &w)); 1928 totpoints = 1; 1929 x[0] = 0.0; 1930 for (PetscInt c = 0; c < Nc; ++c) w[c] = 1.0; 1931 break; 1932 case 1: 1933 ct = DM_POLYTOPE_SEGMENT; 1934 PetscCall(PetscMalloc1(npoints, &ww)); 1935 PetscCall(PetscDTGaussQuadrature(npoints, a, b, x, ww)); 1936 for (PetscInt i = 0; i < npoints; ++i) 1937 for (PetscInt c = 0; c < Nc; ++c) w[i * Nc + c] = ww[i]; 1938 PetscCall(PetscFree(ww)); 1939 break; 1940 case 2: 1941 ct = DM_POLYTOPE_QUADRILATERAL; 1942 PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww)); 1943 PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww)); 1944 for (PetscInt i = 0; i < npoints; ++i) { 1945 for (PetscInt j = 0; j < npoints; ++j) { 1946 x[(i * npoints + j) * dim + 0] = xw[i]; 1947 x[(i * npoints + j) * dim + 1] = xw[j]; 1948 for (PetscInt c = 0; c < Nc; ++c) w[(i * npoints + j) * Nc + c] = ww[i] * ww[j]; 1949 } 1950 } 1951 PetscCall(PetscFree2(xw, ww)); 1952 break; 1953 case 3: 1954 ct = DM_POLYTOPE_HEXAHEDRON; 1955 PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww)); 1956 PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww)); 1957 for (PetscInt i = 0; i < npoints; ++i) { 1958 for (PetscInt j = 0; j < npoints; ++j) { 1959 for (PetscInt k = 0; k < npoints; ++k) { 1960 x[((i * npoints + j) * npoints + k) * dim + 0] = xw[i]; 1961 x[((i * npoints + j) * npoints + k) * dim + 1] = xw[j]; 1962 x[((i * npoints + j) * npoints + k) * dim + 2] = xw[k]; 1963 for (PetscInt c = 0; c < Nc; ++c) w[((i * npoints + j) * npoints + k) * Nc + c] = ww[i] * ww[j] * ww[k]; 1964 } 1965 } 1966 } 1967 PetscCall(PetscFree2(xw, ww)); 1968 break; 1969 default: 1970 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %" PetscInt_FMT, dim); 1971 } 1972 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 1973 PetscCall(PetscQuadratureSetCellType(*q, ct)); 1974 PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1)); 1975 PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w)); 1976 PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "GaussTensor")); 1977 PetscFunctionReturn(PETSC_SUCCESS); 1978 } 1979 1980 /*@ 1981 PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex {cite}`karniadakis2005spectral` 1982 1983 Not Collective 1984 1985 Input Parameters: 1986 + dim - The simplex dimension 1987 . Nc - The number of components 1988 . npoints - The number of points in one dimension 1989 . a - left end of interval (often-1) 1990 - b - right end of interval (often +1) 1991 1992 Output Parameter: 1993 . q - A `PetscQuadrature` object 1994 1995 Level: intermediate 1996 1997 Note: 1998 For `dim` == 1, this is Gauss-Legendre quadrature 1999 2000 .seealso: `PetscDTGaussTensorQuadrature()`, `PetscDTGaussQuadrature()` 2001 @*/ 2002 PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 2003 { 2004 DMPolytopeType ct; 2005 PetscInt totpoints; 2006 PetscReal *p1, *w1; 2007 PetscReal *x, *w; 2008 2009 PetscFunctionBegin; 2010 PetscCheck(!(a != -1.0) && !(b != 1.0), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now"); 2011 switch (dim) { 2012 case 0: 2013 ct = DM_POLYTOPE_POINT; 2014 break; 2015 case 1: 2016 ct = DM_POLYTOPE_SEGMENT; 2017 break; 2018 case 2: 2019 ct = DM_POLYTOPE_TRIANGLE; 2020 break; 2021 case 3: 2022 ct = DM_POLYTOPE_TETRAHEDRON; 2023 break; 2024 default: 2025 ct = DM_POLYTOPE_UNKNOWN; 2026 } 2027 totpoints = 1; 2028 for (PetscInt i = 0; i < dim; ++i) totpoints *= npoints; 2029 PetscCall(PetscMalloc1(totpoints * dim, &x)); 2030 PetscCall(PetscMalloc1(totpoints * Nc, &w)); 2031 PetscCall(PetscMalloc2(npoints, &p1, npoints, &w1)); 2032 for (PetscInt i = 0; i < totpoints * Nc; ++i) w[i] = 1.; 2033 for (PetscInt i = 0, totprev = 1, totrem = totpoints / npoints; i < dim; ++i) { 2034 PetscReal mul; 2035 2036 mul = PetscPowReal(2., -i); 2037 PetscCall(PetscDTGaussJacobiQuadrature(npoints, -1., 1., i, 0.0, p1, w1)); 2038 for (PetscInt pt = 0, l = 0; l < totprev; l++) { 2039 for (PetscInt j = 0; j < npoints; j++) { 2040 for (PetscInt m = 0; m < totrem; m++, pt++) { 2041 for (PetscInt k = 0; k < i; k++) x[pt * dim + k] = (x[pt * dim + k] + 1.) * (1. - p1[j]) * 0.5 - 1.; 2042 x[pt * dim + i] = p1[j]; 2043 for (PetscInt c = 0; c < Nc; c++) w[pt * Nc + c] *= mul * w1[j]; 2044 } 2045 } 2046 } 2047 totprev *= npoints; 2048 totrem /= npoints; 2049 } 2050 PetscCall(PetscFree2(p1, w1)); 2051 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 2052 PetscCall(PetscQuadratureSetCellType(*q, ct)); 2053 PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1)); 2054 PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w)); 2055 PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "StroudConical")); 2056 PetscFunctionReturn(PETSC_SUCCESS); 2057 } 2058 2059 static PetscBool MinSymTriQuadCite = PETSC_FALSE; 2060 const char MinSymTriQuadCitation[] = "@article{WitherdenVincent2015,\n" 2061 " title = {On the identification of symmetric quadrature rules for finite element methods},\n" 2062 " journal = {Computers & Mathematics with Applications},\n" 2063 " volume = {69},\n" 2064 " number = {10},\n" 2065 " pages = {1232-1241},\n" 2066 " year = {2015},\n" 2067 " issn = {0898-1221},\n" 2068 " doi = {10.1016/j.camwa.2015.03.017},\n" 2069 " url = {https://www.sciencedirect.com/science/article/pii/S0898122115001224},\n" 2070 " author = {F.D. Witherden and P.E. Vincent},\n" 2071 "}\n"; 2072 2073 #include "petscdttriquadrules.h" 2074 2075 static PetscBool MinSymTetQuadCite = PETSC_FALSE; 2076 const char MinSymTetQuadCitation[] = "@article{JaskowiecSukumar2021\n" 2077 " author = {Jaskowiec, Jan and Sukumar, N.},\n" 2078 " title = {High-order symmetric cubature rules for tetrahedra and pyramids},\n" 2079 " journal = {International Journal for Numerical Methods in Engineering},\n" 2080 " volume = {122},\n" 2081 " number = {1},\n" 2082 " pages = {148-171},\n" 2083 " doi = {10.1002/nme.6528},\n" 2084 " url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6528},\n" 2085 " eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6528},\n" 2086 " year = {2021}\n" 2087 "}\n"; 2088 2089 #include "petscdttetquadrules.h" 2090 2091 static PetscBool DiagSymTriQuadCite = PETSC_FALSE; 2092 const char DiagSymTriQuadCitation[] = "@article{KongMulderVeldhuizen1999,\n" 2093 " title = {Higher-order triangular and tetrahedral finite elements with mass lumping for solving the wave equation},\n" 2094 " journal = {Journal of Engineering Mathematics},\n" 2095 " volume = {35},\n" 2096 " number = {4},\n" 2097 " pages = {405--426},\n" 2098 " year = {1999},\n" 2099 " doi = {10.1023/A:1004420829610},\n" 2100 " url = {https://link.springer.com/article/10.1023/A:1004420829610},\n" 2101 " author = {MJS Chin-Joe-Kong and Wim A Mulder and Marinus Van Veldhuizen},\n" 2102 "}\n"; 2103 2104 #include "petscdttridiagquadrules.h" 2105 2106 // https://en.wikipedia.org/wiki/Partition_(number_theory) 2107 static PetscErrorCode PetscDTPartitionNumber(PetscInt n, PetscInt *p) 2108 { 2109 // sequence A000041 in the OEIS 2110 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}; 2111 PetscInt tabulated_max = PETSC_STATIC_ARRAY_LENGTH(partition) - 1; 2112 2113 PetscFunctionBegin; 2114 PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Partition number not defined for negative number %" PetscInt_FMT, n); 2115 // not implementing the pentagonal number recurrence, we don't need partition numbers for n that high 2116 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); 2117 *p = partition[n]; 2118 PetscFunctionReturn(PETSC_SUCCESS); 2119 } 2120 2121 /*@ 2122 PetscDTSimplexQuadrature - Create a quadrature rule for a simplex that exactly integrates polynomials up to a given degree. 2123 2124 Not Collective 2125 2126 Input Parameters: 2127 + dim - The spatial dimension of the simplex (1 = segment, 2 = triangle, 3 = tetrahedron) 2128 . degree - The largest polynomial degree that is required to be integrated exactly 2129 - type - `PetscDTSimplexQuadratureType` indicating the type of quadrature rule 2130 2131 Output Parameter: 2132 . quad - A `PetscQuadrature` object for integration over the biunit simplex 2133 (defined by the bounds $x_i >= -1$ and $\sum_i x_i <= 2 - d$) that is exact for 2134 polynomials up to the given degree 2135 2136 Level: intermediate 2137 2138 .seealso: `PetscDTSimplexQuadratureType`, `PetscDTGaussQuadrature()`, `PetscDTStroudCononicalQuadrature()`, `PetscQuadrature` 2139 @*/ 2140 PetscErrorCode PetscDTSimplexQuadrature(PetscInt dim, PetscInt degree, PetscDTSimplexQuadratureType type, PetscQuadrature *quad) 2141 { 2142 PetscDTSimplexQuadratureType orig_type = type; 2143 2144 PetscFunctionBegin; 2145 PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative dimension %" PetscInt_FMT, dim); 2146 PetscCheck(degree >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %" PetscInt_FMT, degree); 2147 if (type == PETSCDTSIMPLEXQUAD_DEFAULT) type = PETSCDTSIMPLEXQUAD_MINSYM; 2148 if (type == PETSCDTSIMPLEXQUAD_CONIC || dim < 2) { 2149 PetscInt points_per_dim = (degree + 2) / 2; // ceil((degree + 1) / 2); 2150 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, points_per_dim, -1, 1, quad)); 2151 } else { 2152 DMPolytopeType ct; 2153 PetscInt n = dim + 1; 2154 PetscInt fact = 1; 2155 PetscInt *part, *perm; 2156 PetscInt p = 0; 2157 PetscInt max_degree; 2158 const PetscInt *nodes_per_type = NULL; 2159 const PetscInt *all_num_full_nodes = NULL; 2160 const PetscReal **weights_list = NULL; 2161 const PetscReal **compact_nodes_list = NULL; 2162 const char *citation = NULL; 2163 PetscBool *cited = NULL; 2164 2165 switch (dim) { 2166 case 0: 2167 ct = DM_POLYTOPE_POINT; 2168 break; 2169 case 1: 2170 ct = DM_POLYTOPE_SEGMENT; 2171 break; 2172 case 2: 2173 ct = DM_POLYTOPE_TRIANGLE; 2174 break; 2175 case 3: 2176 ct = DM_POLYTOPE_TETRAHEDRON; 2177 break; 2178 default: 2179 ct = DM_POLYTOPE_UNKNOWN; 2180 } 2181 if (type == PETSCDTSIMPLEXQUAD_MINSYM) { 2182 switch (dim) { 2183 case 2: 2184 cited = &MinSymTriQuadCite; 2185 citation = MinSymTriQuadCitation; 2186 max_degree = PetscDTWVTriQuad_max_degree; 2187 nodes_per_type = PetscDTWVTriQuad_num_orbits; 2188 all_num_full_nodes = PetscDTWVTriQuad_num_nodes; 2189 weights_list = PetscDTWVTriQuad_weights; 2190 compact_nodes_list = PetscDTWVTriQuad_orbits; 2191 break; 2192 case 3: 2193 cited = &MinSymTetQuadCite; 2194 citation = MinSymTetQuadCitation; 2195 max_degree = PetscDTJSTetQuad_max_degree; 2196 nodes_per_type = PetscDTJSTetQuad_num_orbits; 2197 all_num_full_nodes = PetscDTJSTetQuad_num_nodes; 2198 weights_list = PetscDTJSTetQuad_weights; 2199 compact_nodes_list = PetscDTJSTetQuad_orbits; 2200 break; 2201 default: 2202 max_degree = -1; 2203 break; 2204 } 2205 } else { 2206 switch (dim) { 2207 case 2: 2208 cited = &DiagSymTriQuadCite; 2209 citation = DiagSymTriQuadCitation; 2210 max_degree = PetscDTKMVTriQuad_max_degree; 2211 nodes_per_type = PetscDTKMVTriQuad_num_orbits; 2212 all_num_full_nodes = PetscDTKMVTriQuad_num_nodes; 2213 weights_list = PetscDTKMVTriQuad_weights; 2214 compact_nodes_list = PetscDTKMVTriQuad_orbits; 2215 break; 2216 default: 2217 max_degree = -1; 2218 break; 2219 } 2220 } 2221 2222 if (degree > max_degree) { 2223 PetscCheck(orig_type == PETSCDTSIMPLEXQUAD_DEFAULT, 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); 2224 // fall back to conic 2225 PetscCall(PetscDTSimplexQuadrature(dim, degree, PETSCDTSIMPLEXQUAD_CONIC, quad)); 2226 PetscFunctionReturn(PETSC_SUCCESS); 2227 } 2228 2229 PetscCall(PetscCitationsRegister(citation, cited)); 2230 2231 PetscCall(PetscDTPartitionNumber(n, &p)); 2232 for (PetscInt d = 2; d <= n; d++) fact *= d; 2233 2234 PetscInt num_full_nodes = all_num_full_nodes[degree]; 2235 const PetscReal *all_compact_nodes = compact_nodes_list[degree]; 2236 const PetscReal *all_compact_weights = weights_list[degree]; 2237 nodes_per_type = &nodes_per_type[p * degree]; 2238 2239 PetscReal *points; 2240 PetscReal *counts; 2241 PetscReal *weights; 2242 PetscReal *bary_to_biunit; // row-major transformation of barycentric coordinate to biunit 2243 PetscQuadrature q; 2244 2245 // compute the transformation 2246 PetscCall(PetscMalloc1(n * dim, &bary_to_biunit)); 2247 for (PetscInt d = 0; d < dim; d++) { 2248 for (PetscInt b = 0; b < n; b++) bary_to_biunit[d * n + b] = (d == b) ? 1.0 : -1.0; 2249 } 2250 2251 PetscCall(PetscMalloc3(n, &part, n, &perm, n, &counts)); 2252 PetscCall(PetscCalloc1(num_full_nodes * dim, &points)); 2253 PetscCall(PetscMalloc1(num_full_nodes, &weights)); 2254 2255 // (0, 0, ...) is the first partition lexicographically 2256 PetscCall(PetscArrayzero(part, n)); 2257 PetscCall(PetscArrayzero(counts, n)); 2258 counts[0] = n; 2259 2260 // for each partition 2261 for (PetscInt s = 0, node_offset = 0; s < p; s++) { 2262 PetscInt num_compact_coords = part[n - 1] + 1; 2263 2264 const PetscReal *compact_nodes = all_compact_nodes; 2265 const PetscReal *compact_weights = all_compact_weights; 2266 all_compact_nodes += num_compact_coords * nodes_per_type[s]; 2267 all_compact_weights += nodes_per_type[s]; 2268 2269 // for every permutation of the vertices 2270 for (PetscInt f = 0; f < fact; f++) { 2271 PetscCall(PetscDTEnumPerm(n, f, perm, NULL)); 2272 2273 // check if it is a valid permutation 2274 PetscInt digit; 2275 for (digit = 1; digit < n; digit++) { 2276 // skip permutations that would duplicate a node because it has a smaller symmetry group 2277 if (part[digit - 1] == part[digit] && perm[digit - 1] > perm[digit]) break; 2278 } 2279 if (digit < n) continue; 2280 2281 // create full nodes from this permutation of the compact nodes 2282 PetscReal *full_nodes = &points[node_offset * dim]; 2283 PetscReal *full_weights = &weights[node_offset]; 2284 2285 PetscCall(PetscArraycpy(full_weights, compact_weights, nodes_per_type[s])); 2286 for (PetscInt b = 0; b < n; b++) { 2287 for (PetscInt d = 0; d < dim; d++) { 2288 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]]; 2289 } 2290 } 2291 node_offset += nodes_per_type[s]; 2292 } 2293 2294 if (s < p - 1) { // Generate the next partition 2295 /* A partition is described by the number of coordinates that are in 2296 * each set of duplicates (counts) and redundantly by mapping each 2297 * index to its set of duplicates (part) 2298 * 2299 * Counts should always be in nonincreasing order 2300 * 2301 * We want to generate the partitions lexically by part, which means 2302 * finding the last index where count > 1 and reducing by 1. 2303 * 2304 * For the new counts beyond that index, we eagerly assign the remaining 2305 * capacity of the partition to smaller indices (ensures lexical ordering), 2306 * while respecting the nonincreasing invariant of the counts 2307 */ 2308 PetscInt last_digit = part[n - 1]; 2309 PetscInt last_digit_with_extra = last_digit; 2310 while (counts[last_digit_with_extra] == 1) last_digit_with_extra--; 2311 PetscInt limit = --counts[last_digit_with_extra]; 2312 PetscInt total_to_distribute = last_digit - last_digit_with_extra + 1; 2313 for (PetscInt digit = last_digit_with_extra + 1; digit < n; digit++) { 2314 counts[digit] = PetscMin(limit, total_to_distribute); 2315 total_to_distribute -= PetscMin(limit, total_to_distribute); 2316 } 2317 for (PetscInt digit = 0, offset = 0; digit < n; digit++) { 2318 PetscInt count = counts[digit]; 2319 for (PetscInt c = 0; c < count; c++) part[offset++] = digit; 2320 } 2321 } 2322 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); 2323 } 2324 PetscCall(PetscFree3(part, perm, counts)); 2325 PetscCall(PetscFree(bary_to_biunit)); 2326 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &q)); 2327 PetscCall(PetscQuadratureSetCellType(q, ct)); 2328 PetscCall(PetscQuadratureSetOrder(q, degree)); 2329 PetscCall(PetscQuadratureSetData(q, dim, 1, num_full_nodes, points, weights)); 2330 *quad = q; 2331 } 2332 PetscFunctionReturn(PETSC_SUCCESS); 2333 } 2334 2335 /*@ 2336 PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell 2337 2338 Not Collective 2339 2340 Input Parameters: 2341 + dim - The cell dimension 2342 . level - The number of points in one dimension, $2^l$ 2343 . a - left end of interval (often-1) 2344 - b - right end of interval (often +1) 2345 2346 Output Parameter: 2347 . q - A `PetscQuadrature` object 2348 2349 Level: intermediate 2350 2351 .seealso: `PetscDTGaussTensorQuadrature()`, `PetscQuadrature` 2352 @*/ 2353 PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q) 2354 { 2355 DMPolytopeType ct; 2356 const PetscInt p = 16; /* Digits of precision in the evaluation */ 2357 const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */ 2358 const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */ 2359 const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */ 2360 PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */ 2361 PetscReal wk = 0.5 * PETSC_PI; /* Quadrature weight at x_k */ 2362 PetscReal *x, *w; 2363 PetscInt K, k, npoints; 2364 2365 PetscFunctionBegin; 2366 PetscCheck(dim <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not yet implemented", dim); 2367 PetscCheck(level, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits"); 2368 switch (dim) { 2369 case 0: 2370 ct = DM_POLYTOPE_POINT; 2371 break; 2372 case 1: 2373 ct = DM_POLYTOPE_SEGMENT; 2374 break; 2375 case 2: 2376 ct = DM_POLYTOPE_QUADRILATERAL; 2377 break; 2378 case 3: 2379 ct = DM_POLYTOPE_HEXAHEDRON; 2380 break; 2381 default: 2382 ct = DM_POLYTOPE_UNKNOWN; 2383 } 2384 /* Find K such that the weights are < 32 digits of precision */ 2385 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))); 2386 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 2387 PetscCall(PetscQuadratureSetCellType(*q, ct)); 2388 PetscCall(PetscQuadratureSetOrder(*q, 2 * K + 1)); 2389 npoints = 2 * K - 1; 2390 PetscCall(PetscMalloc1(npoints * dim, &x)); 2391 PetscCall(PetscMalloc1(npoints, &w)); 2392 /* Center term */ 2393 x[0] = beta; 2394 w[0] = 0.5 * alpha * PETSC_PI; 2395 for (k = 1; k < K; ++k) { 2396 wk = 0.5 * alpha * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2397 xk = PetscTanhReal(0.5 * PETSC_PI * PetscSinhReal(k * h)); 2398 x[2 * k - 1] = -alpha * xk + beta; 2399 w[2 * k - 1] = wk; 2400 x[2 * k + 0] = alpha * xk + beta; 2401 w[2 * k + 0] = wk; 2402 } 2403 PetscCall(PetscQuadratureSetData(*q, dim, 1, npoints, x, w)); 2404 PetscFunctionReturn(PETSC_SUCCESS); 2405 } 2406 2407 PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2408 { 2409 const PetscInt p = 16; /* Digits of precision in the evaluation */ 2410 const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */ 2411 const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */ 2412 PetscReal h = 1.0; /* Step size, length between x_k */ 2413 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 2414 PetscReal osum = 0.0; /* Integral on last level */ 2415 PetscReal psum = 0.0; /* Integral on the level before the last level */ 2416 PetscReal sum; /* Integral on current level */ 2417 PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 2418 PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 2419 PetscReal wk; /* Quadrature weight at x_k */ 2420 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 2421 PetscInt d; /* Digits of precision in the integral */ 2422 2423 PetscFunctionBegin; 2424 PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 2425 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 2426 /* Center term */ 2427 func(&beta, ctx, &lval); 2428 sum = 0.5 * alpha * PETSC_PI * lval; 2429 /* */ 2430 do { 2431 PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4; 2432 PetscInt k = 1; 2433 2434 ++l; 2435 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */ 2436 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 2437 psum = osum; 2438 osum = sum; 2439 h *= 0.5; 2440 sum *= 0.5; 2441 do { 2442 wk = 0.5 * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2443 yk = 1.0 / (PetscExpReal(0.5 * PETSC_PI * PetscSinhReal(k * h)) * PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2444 lx = -alpha * (1.0 - yk) + beta; 2445 rx = alpha * (1.0 - yk) + beta; 2446 func(&lx, ctx, &lval); 2447 func(&rx, ctx, &rval); 2448 lterm = alpha * wk * lval; 2449 maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm); 2450 sum += lterm; 2451 rterm = alpha * wk * rval; 2452 maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm); 2453 sum += rterm; 2454 ++k; 2455 /* Only need to evaluate every other point on refined levels */ 2456 if (l != 1) ++k; 2457 } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */ 2458 2459 d1 = PetscLog10Real(PetscAbsReal(sum - osum)); 2460 d2 = PetscLog10Real(PetscAbsReal(sum - psum)); 2461 d3 = PetscLog10Real(maxTerm) - p; 2462 if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0; 2463 else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm))); 2464 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4))); 2465 } while (d < digits && l < 12); 2466 *sol = sum; 2467 PetscCall(PetscFPTrapPop()); 2468 PetscFunctionReturn(PETSC_SUCCESS); 2469 } 2470 2471 #if defined(PETSC_HAVE_MPFR) 2472 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2473 { 2474 const PetscInt safetyFactor = 2; /* Calculate abscissa until 2*p digits */ 2475 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 2476 mpfr_t alpha; /* Half-width of the integration interval */ 2477 mpfr_t beta; /* Center of the integration interval */ 2478 mpfr_t h; /* Step size, length between x_k */ 2479 mpfr_t osum; /* Integral on last level */ 2480 mpfr_t psum; /* Integral on the level before the last level */ 2481 mpfr_t sum; /* Integral on current level */ 2482 mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 2483 mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 2484 mpfr_t wk; /* Quadrature weight at x_k */ 2485 PetscReal lval, rval, rtmp; /* Terms in the quadature sum to the left and right of 0 */ 2486 PetscInt d; /* Digits of precision in the integral */ 2487 mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp; 2488 2489 PetscFunctionBegin; 2490 PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 2491 /* Create high precision storage */ 2492 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); 2493 /* Initialization */ 2494 mpfr_set_d(alpha, 0.5 * (b - a), MPFR_RNDN); 2495 mpfr_set_d(beta, 0.5 * (b + a), MPFR_RNDN); 2496 mpfr_set_d(osum, 0.0, MPFR_RNDN); 2497 mpfr_set_d(psum, 0.0, MPFR_RNDN); 2498 mpfr_set_d(h, 1.0, MPFR_RNDN); 2499 mpfr_const_pi(pi2, MPFR_RNDN); 2500 mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN); 2501 /* Center term */ 2502 rtmp = 0.5 * (b + a); 2503 func(&rtmp, ctx, &lval); 2504 mpfr_set(sum, pi2, MPFR_RNDN); 2505 mpfr_mul(sum, sum, alpha, MPFR_RNDN); 2506 mpfr_mul_d(sum, sum, lval, MPFR_RNDN); 2507 /* */ 2508 do { 2509 PetscReal d1, d2, d3, d4; 2510 PetscInt k = 1; 2511 2512 ++l; 2513 mpfr_set_d(maxTerm, 0.0, MPFR_RNDN); 2514 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */ 2515 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 2516 mpfr_set(psum, osum, MPFR_RNDN); 2517 mpfr_set(osum, sum, MPFR_RNDN); 2518 mpfr_mul_d(h, h, 0.5, MPFR_RNDN); 2519 mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN); 2520 do { 2521 mpfr_set_si(kh, k, MPFR_RNDN); 2522 mpfr_mul(kh, kh, h, MPFR_RNDN); 2523 /* Weight */ 2524 mpfr_set(wk, h, MPFR_RNDN); 2525 mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN); 2526 mpfr_mul(msinh, msinh, pi2, MPFR_RNDN); 2527 mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN); 2528 mpfr_cosh(tmp, msinh, MPFR_RNDN); 2529 mpfr_sqr(tmp, tmp, MPFR_RNDN); 2530 mpfr_mul(wk, wk, mcosh, MPFR_RNDN); 2531 mpfr_div(wk, wk, tmp, MPFR_RNDN); 2532 /* Abscissa */ 2533 mpfr_set_d(yk, 1.0, MPFR_RNDZ); 2534 mpfr_cosh(tmp, msinh, MPFR_RNDN); 2535 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 2536 mpfr_exp(tmp, msinh, MPFR_RNDN); 2537 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 2538 /* Quadrature points */ 2539 mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ); 2540 mpfr_mul(lx, lx, alpha, MPFR_RNDU); 2541 mpfr_add(lx, lx, beta, MPFR_RNDU); 2542 mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ); 2543 mpfr_mul(rx, rx, alpha, MPFR_RNDD); 2544 mpfr_add(rx, rx, beta, MPFR_RNDD); 2545 /* Evaluation */ 2546 rtmp = mpfr_get_d(lx, MPFR_RNDU); 2547 func(&rtmp, ctx, &lval); 2548 rtmp = mpfr_get_d(rx, MPFR_RNDD); 2549 func(&rtmp, ctx, &rval); 2550 /* Update */ 2551 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 2552 mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN); 2553 mpfr_add(sum, sum, tmp, MPFR_RNDN); 2554 mpfr_abs(tmp, tmp, MPFR_RNDN); 2555 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 2556 mpfr_set(curTerm, tmp, MPFR_RNDN); 2557 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 2558 mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN); 2559 mpfr_add(sum, sum, tmp, MPFR_RNDN); 2560 mpfr_abs(tmp, tmp, MPFR_RNDN); 2561 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 2562 mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN); 2563 ++k; 2564 /* Only need to evaluate every other point on refined levels */ 2565 if (l != 1) ++k; 2566 mpfr_log10(tmp, wk, MPFR_RNDN); 2567 mpfr_abs(tmp, tmp, MPFR_RNDN); 2568 } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor * digits); /* Only need to evaluate sum until weights are < 32 digits of precision */ 2569 mpfr_sub(tmp, sum, osum, MPFR_RNDN); 2570 mpfr_abs(tmp, tmp, MPFR_RNDN); 2571 mpfr_log10(tmp, tmp, MPFR_RNDN); 2572 d1 = mpfr_get_d(tmp, MPFR_RNDN); 2573 mpfr_sub(tmp, sum, psum, MPFR_RNDN); 2574 mpfr_abs(tmp, tmp, MPFR_RNDN); 2575 mpfr_log10(tmp, tmp, MPFR_RNDN); 2576 d2 = mpfr_get_d(tmp, MPFR_RNDN); 2577 mpfr_log10(tmp, maxTerm, MPFR_RNDN); 2578 d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits; 2579 mpfr_log10(tmp, curTerm, MPFR_RNDN); 2580 d4 = mpfr_get_d(tmp, MPFR_RNDN); 2581 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4))); 2582 } while (d < digits && l < 8); 2583 *sol = mpfr_get_d(sum, MPFR_RNDN); 2584 /* Cleanup */ 2585 mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL); 2586 PetscFunctionReturn(PETSC_SUCCESS); 2587 } 2588 #else 2589 2590 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2591 { 2592 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp"); 2593 } 2594 #endif 2595 2596 /*@ 2597 PetscDTTensorQuadratureCreate - create the tensor product quadrature from two lower-dimensional quadratures 2598 2599 Not Collective 2600 2601 Input Parameters: 2602 + q1 - The first quadrature 2603 - q2 - The second quadrature 2604 2605 Output Parameter: 2606 . q - A `PetscQuadrature` object 2607 2608 Level: intermediate 2609 2610 .seealso: `PetscQuadrature`, `PetscDTGaussTensorQuadrature()` 2611 @*/ 2612 PetscErrorCode PetscDTTensorQuadratureCreate(PetscQuadrature q1, PetscQuadrature q2, PetscQuadrature *q) 2613 { 2614 DMPolytopeType ct1, ct2, ct; 2615 const PetscReal *x1, *w1, *x2, *w2; 2616 PetscReal *x, *w; 2617 PetscInt dim1, Nc1, Np1, order1, qa, d1; 2618 PetscInt dim2, Nc2, Np2, order2, qb, d2; 2619 PetscInt dim, Nc, Np, order, qc, d; 2620 2621 PetscFunctionBegin; 2622 PetscValidHeaderSpecific(q1, PETSCQUADRATURE_CLASSID, 1); 2623 PetscValidHeaderSpecific(q2, PETSCQUADRATURE_CLASSID, 2); 2624 PetscAssertPointer(q, 3); 2625 2626 PetscCall(PetscQuadratureGetOrder(q1, &order1)); 2627 PetscCall(PetscQuadratureGetOrder(q2, &order2)); 2628 PetscCheck(order1 == order2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Order1 %" PetscInt_FMT " != %" PetscInt_FMT " Order2", order1, order2); 2629 PetscCall(PetscQuadratureGetData(q1, &dim1, &Nc1, &Np1, &x1, &w1)); 2630 PetscCall(PetscQuadratureGetCellType(q1, &ct1)); 2631 PetscCall(PetscQuadratureGetData(q2, &dim2, &Nc2, &Np2, &x2, &w2)); 2632 PetscCall(PetscQuadratureGetCellType(q2, &ct2)); 2633 PetscCheck(Nc1 == Nc2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "NumComp1 %" PetscInt_FMT " != %" PetscInt_FMT " NumComp2", Nc1, Nc2); 2634 2635 switch (ct1) { 2636 case DM_POLYTOPE_POINT: 2637 ct = ct2; 2638 break; 2639 case DM_POLYTOPE_SEGMENT: 2640 switch (ct2) { 2641 case DM_POLYTOPE_POINT: 2642 ct = ct1; 2643 break; 2644 case DM_POLYTOPE_SEGMENT: 2645 ct = DM_POLYTOPE_QUADRILATERAL; 2646 break; 2647 case DM_POLYTOPE_TRIANGLE: 2648 ct = DM_POLYTOPE_TRI_PRISM; 2649 break; 2650 case DM_POLYTOPE_QUADRILATERAL: 2651 ct = DM_POLYTOPE_HEXAHEDRON; 2652 break; 2653 case DM_POLYTOPE_TETRAHEDRON: 2654 ct = DM_POLYTOPE_UNKNOWN; 2655 break; 2656 case DM_POLYTOPE_HEXAHEDRON: 2657 ct = DM_POLYTOPE_UNKNOWN; 2658 break; 2659 default: 2660 ct = DM_POLYTOPE_UNKNOWN; 2661 } 2662 break; 2663 case DM_POLYTOPE_TRIANGLE: 2664 switch (ct2) { 2665 case DM_POLYTOPE_POINT: 2666 ct = ct1; 2667 break; 2668 case DM_POLYTOPE_SEGMENT: 2669 ct = DM_POLYTOPE_TRI_PRISM; 2670 break; 2671 case DM_POLYTOPE_TRIANGLE: 2672 ct = DM_POLYTOPE_UNKNOWN; 2673 break; 2674 case DM_POLYTOPE_QUADRILATERAL: 2675 ct = DM_POLYTOPE_UNKNOWN; 2676 break; 2677 case DM_POLYTOPE_TETRAHEDRON: 2678 ct = DM_POLYTOPE_UNKNOWN; 2679 break; 2680 case DM_POLYTOPE_HEXAHEDRON: 2681 ct = DM_POLYTOPE_UNKNOWN; 2682 break; 2683 default: 2684 ct = DM_POLYTOPE_UNKNOWN; 2685 } 2686 break; 2687 case DM_POLYTOPE_QUADRILATERAL: 2688 switch (ct2) { 2689 case DM_POLYTOPE_POINT: 2690 ct = ct1; 2691 break; 2692 case DM_POLYTOPE_SEGMENT: 2693 ct = DM_POLYTOPE_HEXAHEDRON; 2694 break; 2695 case DM_POLYTOPE_TRIANGLE: 2696 ct = DM_POLYTOPE_UNKNOWN; 2697 break; 2698 case DM_POLYTOPE_QUADRILATERAL: 2699 ct = DM_POLYTOPE_UNKNOWN; 2700 break; 2701 case DM_POLYTOPE_TETRAHEDRON: 2702 ct = DM_POLYTOPE_UNKNOWN; 2703 break; 2704 case DM_POLYTOPE_HEXAHEDRON: 2705 ct = DM_POLYTOPE_UNKNOWN; 2706 break; 2707 default: 2708 ct = DM_POLYTOPE_UNKNOWN; 2709 } 2710 break; 2711 case DM_POLYTOPE_TETRAHEDRON: 2712 switch (ct2) { 2713 case DM_POLYTOPE_POINT: 2714 ct = ct1; 2715 break; 2716 case DM_POLYTOPE_SEGMENT: 2717 ct = DM_POLYTOPE_UNKNOWN; 2718 break; 2719 case DM_POLYTOPE_TRIANGLE: 2720 ct = DM_POLYTOPE_UNKNOWN; 2721 break; 2722 case DM_POLYTOPE_QUADRILATERAL: 2723 ct = DM_POLYTOPE_UNKNOWN; 2724 break; 2725 case DM_POLYTOPE_TETRAHEDRON: 2726 ct = DM_POLYTOPE_UNKNOWN; 2727 break; 2728 case DM_POLYTOPE_HEXAHEDRON: 2729 ct = DM_POLYTOPE_UNKNOWN; 2730 break; 2731 default: 2732 ct = DM_POLYTOPE_UNKNOWN; 2733 } 2734 break; 2735 case DM_POLYTOPE_HEXAHEDRON: 2736 switch (ct2) { 2737 case DM_POLYTOPE_POINT: 2738 ct = ct1; 2739 break; 2740 case DM_POLYTOPE_SEGMENT: 2741 ct = DM_POLYTOPE_UNKNOWN; 2742 break; 2743 case DM_POLYTOPE_TRIANGLE: 2744 ct = DM_POLYTOPE_UNKNOWN; 2745 break; 2746 case DM_POLYTOPE_QUADRILATERAL: 2747 ct = DM_POLYTOPE_UNKNOWN; 2748 break; 2749 case DM_POLYTOPE_TETRAHEDRON: 2750 ct = DM_POLYTOPE_UNKNOWN; 2751 break; 2752 case DM_POLYTOPE_HEXAHEDRON: 2753 ct = DM_POLYTOPE_UNKNOWN; 2754 break; 2755 default: 2756 ct = DM_POLYTOPE_UNKNOWN; 2757 } 2758 break; 2759 default: 2760 ct = DM_POLYTOPE_UNKNOWN; 2761 } 2762 dim = dim1 + dim2; 2763 Nc = Nc1; 2764 Np = Np1 * Np2; 2765 order = order1; 2766 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 2767 PetscCall(PetscQuadratureSetCellType(*q, ct)); 2768 PetscCall(PetscQuadratureSetOrder(*q, order)); 2769 PetscCall(PetscMalloc1(Np * dim, &x)); 2770 PetscCall(PetscMalloc1(Np, &w)); 2771 for (qa = 0, qc = 0; qa < Np1; ++qa) { 2772 for (qb = 0; qb < Np2; ++qb, ++qc) { 2773 for (d1 = 0, d = 0; d1 < dim1; ++d1, ++d) x[qc * dim + d] = x1[qa * dim1 + d1]; 2774 for (d2 = 0; d2 < dim2; ++d2, ++d) x[qc * dim + d] = x2[qb * dim2 + d2]; 2775 w[qc] = w1[qa] * w2[qb]; 2776 } 2777 } 2778 PetscCall(PetscQuadratureSetData(*q, dim, Nc, Np, x, w)); 2779 PetscFunctionReturn(PETSC_SUCCESS); 2780 } 2781 2782 /* Overwrites A. Can only handle full-rank problems with m>=n 2783 A in column-major format 2784 Ainv in row-major format 2785 tau has length m 2786 worksize must be >= max(1,n) 2787 */ 2788 static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m, PetscInt mstride, PetscInt n, PetscReal *A_in, PetscReal *Ainv_out, PetscScalar *tau, PetscInt worksize, PetscScalar *work) 2789 { 2790 PetscBLASInt M, N, K, lda, ldb, ldwork, info; 2791 PetscScalar *A, *Ainv, *R, *Q, Alpha; 2792 2793 PetscFunctionBegin; 2794 #if defined(PETSC_USE_COMPLEX) 2795 { 2796 PetscInt i, j; 2797 PetscCall(PetscMalloc2(m * n, &A, m * n, &Ainv)); 2798 for (j = 0; j < n; j++) { 2799 for (i = 0; i < m; i++) A[i + m * j] = A_in[i + mstride * j]; 2800 } 2801 mstride = m; 2802 } 2803 #else 2804 A = A_in; 2805 Ainv = Ainv_out; 2806 #endif 2807 2808 PetscCall(PetscBLASIntCast(m, &M)); 2809 PetscCall(PetscBLASIntCast(n, &N)); 2810 PetscCall(PetscBLASIntCast(mstride, &lda)); 2811 PetscCall(PetscBLASIntCast(worksize, &ldwork)); 2812 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 2813 PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info)); 2814 PetscCall(PetscFPTrapPop()); 2815 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error"); 2816 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */ 2817 2818 /* Extract an explicit representation of Q */ 2819 Q = Ainv; 2820 PetscCall(PetscArraycpy(Q, A, mstride * n)); 2821 K = N; /* full rank */ 2822 PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info)); 2823 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error"); 2824 2825 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */ 2826 Alpha = 1.0; 2827 ldb = lda; 2828 PetscCallBLAS("BLAStrsm", BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb)); 2829 /* Ainv is Q, overwritten with inverse */ 2830 2831 #if defined(PETSC_USE_COMPLEX) 2832 { 2833 PetscInt i; 2834 for (i = 0; i < m * n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]); 2835 PetscCall(PetscFree2(A, Ainv)); 2836 } 2837 #endif 2838 PetscFunctionReturn(PETSC_SUCCESS); 2839 } 2840 2841 /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */ 2842 static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval, const PetscReal *x, PetscInt ndegree, const PetscInt *degrees, PetscBool Transpose, PetscReal *B) 2843 { 2844 PetscReal *Bv; 2845 PetscInt i, j; 2846 2847 PetscFunctionBegin; 2848 PetscCall(PetscMalloc1((ninterval + 1) * ndegree, &Bv)); 2849 /* Point evaluation of L_p on all the source vertices */ 2850 PetscCall(PetscDTLegendreEval(ninterval + 1, x, ndegree, degrees, Bv, NULL, NULL)); 2851 /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */ 2852 for (i = 0; i < ninterval; i++) { 2853 for (j = 0; j < ndegree; j++) { 2854 if (Transpose) B[i + ninterval * j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j]; 2855 else B[i * ndegree + j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j]; 2856 } 2857 } 2858 PetscCall(PetscFree(Bv)); 2859 PetscFunctionReturn(PETSC_SUCCESS); 2860 } 2861 2862 /*@ 2863 PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals 2864 2865 Not Collective 2866 2867 Input Parameters: 2868 + degree - degree of reconstruction polynomial 2869 . nsource - number of source intervals 2870 . sourcex - sorted coordinates of source cell boundaries (length `nsource`+1) 2871 . ntarget - number of target intervals 2872 - targetx - sorted coordinates of target cell boundaries (length `ntarget`+1) 2873 2874 Output Parameter: 2875 . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s] 2876 2877 Level: advanced 2878 2879 .seealso: `PetscDTLegendreEval()` 2880 @*/ 2881 PetscErrorCode PetscDTReconstructPoly(PetscInt degree, PetscInt nsource, const PetscReal sourcex[], PetscInt ntarget, const PetscReal targetx[], PetscReal R[]) 2882 { 2883 PetscInt i, j, k, *bdegrees, worksize; 2884 PetscReal xmin, xmax, center, hscale, *sourcey, *targety, *Bsource, *Bsinv, *Btarget; 2885 PetscScalar *tau, *work; 2886 2887 PetscFunctionBegin; 2888 PetscAssertPointer(sourcex, 3); 2889 PetscAssertPointer(targetx, 5); 2890 PetscAssertPointer(R, 6); 2891 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); 2892 if (PetscDefined(USE_DEBUG)) { 2893 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]); 2894 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]); 2895 } 2896 xmin = PetscMin(sourcex[0], targetx[0]); 2897 xmax = PetscMax(sourcex[nsource], targetx[ntarget]); 2898 center = (xmin + xmax) / 2; 2899 hscale = (xmax - xmin) / 2; 2900 worksize = nsource; 2901 PetscCall(PetscMalloc4(degree + 1, &bdegrees, nsource + 1, &sourcey, nsource * (degree + 1), &Bsource, worksize, &work)); 2902 PetscCall(PetscMalloc4(nsource, &tau, nsource * (degree + 1), &Bsinv, ntarget + 1, &targety, ntarget * (degree + 1), &Btarget)); 2903 for (i = 0; i <= nsource; i++) sourcey[i] = (sourcex[i] - center) / hscale; 2904 for (i = 0; i <= degree; i++) bdegrees[i] = i + 1; 2905 PetscCall(PetscDTLegendreIntegrate(nsource, sourcey, degree + 1, bdegrees, PETSC_TRUE, Bsource)); 2906 PetscCall(PetscDTPseudoInverseQR(nsource, nsource, degree + 1, Bsource, Bsinv, tau, nsource, work)); 2907 for (i = 0; i <= ntarget; i++) targety[i] = (targetx[i] - center) / hscale; 2908 PetscCall(PetscDTLegendreIntegrate(ntarget, targety, degree + 1, bdegrees, PETSC_FALSE, Btarget)); 2909 for (i = 0; i < ntarget; i++) { 2910 PetscReal rowsum = 0; 2911 for (j = 0; j < nsource; j++) { 2912 PetscReal sum = 0; 2913 for (k = 0; k < degree + 1; k++) sum += Btarget[i * (degree + 1) + k] * Bsinv[k * nsource + j]; 2914 R[i * nsource + j] = sum; 2915 rowsum += sum; 2916 } 2917 for (j = 0; j < nsource; j++) R[i * nsource + j] /= rowsum; /* normalize each row */ 2918 } 2919 PetscCall(PetscFree4(bdegrees, sourcey, Bsource, work)); 2920 PetscCall(PetscFree4(tau, Bsinv, targety, Btarget)); 2921 PetscFunctionReturn(PETSC_SUCCESS); 2922 } 2923 2924 /*@ 2925 PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points 2926 2927 Not Collective 2928 2929 Input Parameters: 2930 + n - the number of GLL nodes 2931 . nodes - the GLL nodes 2932 . weights - the GLL weights 2933 - f - the function values at the nodes 2934 2935 Output Parameter: 2936 . in - the value of the integral 2937 2938 Level: beginner 2939 2940 .seealso: `PetscDTGaussLobattoLegendreQuadrature()` 2941 @*/ 2942 PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n, PetscReal nodes[], PetscReal weights[], const PetscReal f[], PetscReal *in) 2943 { 2944 PetscInt i; 2945 2946 PetscFunctionBegin; 2947 *in = 0.; 2948 for (i = 0; i < n; i++) *in += f[i] * f[i] * weights[i]; 2949 PetscFunctionReturn(PETSC_SUCCESS); 2950 } 2951 2952 /*@C 2953 PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element 2954 2955 Not Collective 2956 2957 Input Parameters: 2958 + n - the number of GLL nodes 2959 . nodes - the GLL nodes, of length `n` 2960 - weights - the GLL weights, of length `n` 2961 2962 Output Parameter: 2963 . AA - the stiffness element, of size `n` by `n` 2964 2965 Level: beginner 2966 2967 Notes: 2968 Destroy this with `PetscGaussLobattoLegendreElementLaplacianDestroy()` 2969 2970 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) 2971 2972 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()` 2973 @*/ 2974 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 2975 { 2976 PetscReal **A; 2977 const PetscReal *gllnodes = nodes; 2978 const PetscInt p = n - 1; 2979 PetscReal z0, z1, z2 = -1, x, Lpj, Lpr; 2980 PetscInt i, j, nn, r; 2981 2982 PetscFunctionBegin; 2983 PetscCall(PetscMalloc1(n, &A)); 2984 PetscCall(PetscMalloc1(n * n, &A[0])); 2985 for (i = 1; i < n; i++) A[i] = A[i - 1] + n; 2986 2987 for (j = 1; j < p; j++) { 2988 x = gllnodes[j]; 2989 z0 = 1.; 2990 z1 = x; 2991 for (nn = 1; nn < p; nn++) { 2992 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 2993 z0 = z1; 2994 z1 = z2; 2995 } 2996 Lpj = z2; 2997 for (r = 1; r < p; r++) { 2998 if (r == j) { 2999 A[j][j] = 2. / (3. * (1. - gllnodes[j] * gllnodes[j]) * Lpj * Lpj); 3000 } else { 3001 x = gllnodes[r]; 3002 z0 = 1.; 3003 z1 = x; 3004 for (nn = 1; nn < p; nn++) { 3005 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 3006 z0 = z1; 3007 z1 = z2; 3008 } 3009 Lpr = z2; 3010 A[r][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * Lpr * (gllnodes[j] - gllnodes[r]) * (gllnodes[j] - gllnodes[r])); 3011 } 3012 } 3013 } 3014 for (j = 1; j < p + 1; j++) { 3015 x = gllnodes[j]; 3016 z0 = 1.; 3017 z1 = x; 3018 for (nn = 1; nn < p; nn++) { 3019 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 3020 z0 = z1; 3021 z1 = z2; 3022 } 3023 Lpj = z2; 3024 A[j][0] = 4. * PetscPowRealInt(-1., p) / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. + gllnodes[j]) * (1. + gllnodes[j])); 3025 A[0][j] = A[j][0]; 3026 } 3027 for (j = 0; j < p; j++) { 3028 x = gllnodes[j]; 3029 z0 = 1.; 3030 z1 = x; 3031 for (nn = 1; nn < p; nn++) { 3032 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 3033 z0 = z1; 3034 z1 = z2; 3035 } 3036 Lpj = z2; 3037 3038 A[p][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. - gllnodes[j]) * (1. - gllnodes[j])); 3039 A[j][p] = A[p][j]; 3040 } 3041 A[0][0] = 0.5 + (((PetscReal)p) * (((PetscReal)p) + 1.) - 2.) / 6.; 3042 A[p][p] = A[0][0]; 3043 *AA = A; 3044 PetscFunctionReturn(PETSC_SUCCESS); 3045 } 3046 3047 /*@C 3048 PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element created with `PetscGaussLobattoLegendreElementLaplacianCreate()` 3049 3050 Not Collective 3051 3052 Input Parameters: 3053 + n - the number of GLL nodes 3054 . nodes - the GLL nodes, ignored 3055 . weights - the GLL weightss, ignored 3056 - AA - the stiffness element from `PetscGaussLobattoLegendreElementLaplacianCreate()` 3057 3058 Level: beginner 3059 3060 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()` 3061 @*/ 3062 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 3063 { 3064 PetscFunctionBegin; 3065 PetscCall(PetscFree((*AA)[0])); 3066 PetscCall(PetscFree(*AA)); 3067 *AA = NULL; 3068 PetscFunctionReturn(PETSC_SUCCESS); 3069 } 3070 3071 /*@C 3072 PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element 3073 3074 Not Collective 3075 3076 Input Parameters: 3077 + n - the number of GLL nodes 3078 . nodes - the GLL nodes, of length `n` 3079 - weights - the GLL weights, of length `n` 3080 3081 Output Parameters: 3082 + AA - the stiffness element, of dimension `n` by `n` 3083 - AAT - the transpose of AA (pass in `NULL` if you do not need this array), of dimension `n` by `n` 3084 3085 Level: beginner 3086 3087 Notes: 3088 Destroy this with `PetscGaussLobattoLegendreElementGradientDestroy()` 3089 3090 You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row-oriented 3091 3092 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`, `PetscGaussLobattoLegendreElementGradientDestroy()` 3093 @*/ 3094 PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA, PetscReal ***AAT) 3095 { 3096 PetscReal **A, **AT = NULL; 3097 const PetscReal *gllnodes = nodes; 3098 const PetscInt p = n - 1; 3099 PetscReal Li, Lj, d0; 3100 PetscInt i, j; 3101 3102 PetscFunctionBegin; 3103 PetscCall(PetscMalloc1(n, &A)); 3104 PetscCall(PetscMalloc1(n * n, &A[0])); 3105 for (i = 1; i < n; i++) A[i] = A[i - 1] + n; 3106 3107 if (AAT) { 3108 PetscCall(PetscMalloc1(n, &AT)); 3109 PetscCall(PetscMalloc1(n * n, &AT[0])); 3110 for (i = 1; i < n; i++) AT[i] = AT[i - 1] + n; 3111 } 3112 3113 if (n == 1) A[0][0] = 0.; 3114 d0 = (PetscReal)p * ((PetscReal)p + 1.) / 4.; 3115 for (i = 0; i < n; i++) { 3116 for (j = 0; j < n; j++) { 3117 A[i][j] = 0.; 3118 PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li)); 3119 PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj)); 3120 if (i != j) A[i][j] = Li / (Lj * (gllnodes[i] - gllnodes[j])); 3121 if ((j == i) && (i == 0)) A[i][j] = -d0; 3122 if (j == i && i == p) A[i][j] = d0; 3123 if (AT) AT[j][i] = A[i][j]; 3124 } 3125 } 3126 if (AAT) *AAT = AT; 3127 *AA = A; 3128 PetscFunctionReturn(PETSC_SUCCESS); 3129 } 3130 3131 /*@C 3132 PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with `PetscGaussLobattoLegendreElementGradientCreate()` 3133 3134 Not Collective 3135 3136 Input Parameters: 3137 + n - the number of GLL nodes 3138 . nodes - the GLL nodes, ignored 3139 . weights - the GLL weights, ignored 3140 . AA - the stiffness element obtained with `PetscGaussLobattoLegendreElementGradientCreate()` 3141 - AAT - the transpose of the element obtained with `PetscGaussLobattoLegendreElementGradientCreate()` 3142 3143 Level: beginner 3144 3145 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionCreate()` 3146 @*/ 3147 PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA, PetscReal ***AAT) 3148 { 3149 PetscFunctionBegin; 3150 PetscCall(PetscFree((*AA)[0])); 3151 PetscCall(PetscFree(*AA)); 3152 *AA = NULL; 3153 if (AAT) { 3154 PetscCall(PetscFree((*AAT)[0])); 3155 PetscCall(PetscFree(*AAT)); 3156 *AAT = NULL; 3157 } 3158 PetscFunctionReturn(PETSC_SUCCESS); 3159 } 3160 3161 /*@C 3162 PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element 3163 3164 Not Collective 3165 3166 Input Parameters: 3167 + n - the number of GLL nodes 3168 . nodes - the GLL nodes, of length `n` 3169 - weights - the GLL weights, of length `n` 3170 3171 Output Parameter: 3172 . AA - the stiffness element, of dimension `n` by `n` 3173 3174 Level: beginner 3175 3176 Notes: 3177 Destroy this with `PetscGaussLobattoLegendreElementAdvectionDestroy()` 3178 3179 This is the same as the Gradient operator multiplied by the diagonal mass matrix 3180 3181 You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row-oriented 3182 3183 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionDestroy()` 3184 @*/ 3185 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 3186 { 3187 PetscReal **D; 3188 const PetscReal *gllweights = weights; 3189 const PetscInt glln = n; 3190 PetscInt i, j; 3191 3192 PetscFunctionBegin; 3193 PetscCall(PetscGaussLobattoLegendreElementGradientCreate(n, nodes, weights, &D, NULL)); 3194 for (i = 0; i < glln; i++) { 3195 for (j = 0; j < glln; j++) D[i][j] = gllweights[i] * D[i][j]; 3196 } 3197 *AA = D; 3198 PetscFunctionReturn(PETSC_SUCCESS); 3199 } 3200 3201 /*@C 3202 PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element created with `PetscGaussLobattoLegendreElementAdvectionCreate()` 3203 3204 Not Collective 3205 3206 Input Parameters: 3207 + n - the number of GLL nodes 3208 . nodes - the GLL nodes, ignored 3209 . weights - the GLL weights, ignored 3210 - AA - advection obtained with `PetscGaussLobattoLegendreElementAdvectionCreate()` 3211 3212 Level: beginner 3213 3214 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementAdvectionCreate()` 3215 @*/ 3216 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA) 3217 { 3218 PetscFunctionBegin; 3219 PetscCall(PetscFree((*AA)[0])); 3220 PetscCall(PetscFree(*AA)); 3221 *AA = NULL; 3222 PetscFunctionReturn(PETSC_SUCCESS); 3223 } 3224 3225 PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 3226 { 3227 PetscReal **A; 3228 const PetscReal *gllweights = weights; 3229 const PetscInt glln = n; 3230 PetscInt i, j; 3231 3232 PetscFunctionBegin; 3233 PetscCall(PetscMalloc1(glln, &A)); 3234 PetscCall(PetscMalloc1(glln * glln, &A[0])); 3235 for (i = 1; i < glln; i++) A[i] = A[i - 1] + glln; 3236 if (glln == 1) A[0][0] = 0.; 3237 for (i = 0; i < glln; i++) { 3238 for (j = 0; j < glln; j++) { 3239 A[i][j] = 0.; 3240 if (j == i) A[i][j] = gllweights[i]; 3241 } 3242 } 3243 *AA = A; 3244 PetscFunctionReturn(PETSC_SUCCESS); 3245 } 3246 3247 PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 3248 { 3249 PetscFunctionBegin; 3250 PetscCall(PetscFree((*AA)[0])); 3251 PetscCall(PetscFree(*AA)); 3252 *AA = NULL; 3253 PetscFunctionReturn(PETSC_SUCCESS); 3254 } 3255 3256 /*@ 3257 PetscDTIndexToBary - convert an index into a barycentric coordinate. 3258 3259 Input Parameters: 3260 + 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) 3261 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 3262 - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum) 3263 3264 Output Parameter: 3265 . coord - will be filled with the barycentric coordinate, of length `n` 3266 3267 Level: beginner 3268 3269 Note: 3270 The indices map to barycentric coordinates in lexicographic order, where the first index is the 3271 least significant and the last index is the most significant. 3272 3273 .seealso: `PetscDTBaryToIndex()` 3274 @*/ 3275 PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[]) 3276 { 3277 PetscInt c, d, s, total, subtotal, nexttotal; 3278 3279 PetscFunctionBeginHot; 3280 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 3281 PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative"); 3282 if (!len) { 3283 PetscCheck(!sum && !index, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 3284 PetscFunctionReturn(PETSC_SUCCESS); 3285 } 3286 for (c = 1, total = 1; c <= len; c++) { 3287 /* total is the number of ways to have a tuple of length c with sum */ 3288 if (index < total) break; 3289 total = (total * (sum + c)) / c; 3290 } 3291 PetscCheck(c <= len, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range"); 3292 for (d = c; d < len; d++) coord[d] = 0; 3293 for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) { 3294 /* subtotal is the number of ways to have a tuple of length c with sum s */ 3295 /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */ 3296 if ((index + subtotal) >= total) { 3297 coord[--c] = sum - s; 3298 index -= (total - subtotal); 3299 sum = s; 3300 total = nexttotal; 3301 subtotal = 1; 3302 nexttotal = 1; 3303 s = 0; 3304 } else { 3305 subtotal = (subtotal * (c + s)) / (s + 1); 3306 nexttotal = (nexttotal * (c - 1 + s)) / (s + 1); 3307 s++; 3308 } 3309 } 3310 PetscFunctionReturn(PETSC_SUCCESS); 3311 } 3312 3313 /*@ 3314 PetscDTBaryToIndex - convert a barycentric coordinate to an index 3315 3316 Input Parameters: 3317 + 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) 3318 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 3319 - coord - a barycentric coordinate with the given length `len` and `sum` 3320 3321 Output Parameter: 3322 . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum) 3323 3324 Level: beginner 3325 3326 Note: 3327 The indices map to barycentric coordinates in lexicographic order, where the first index is the 3328 least significant and the last index is the most significant. 3329 3330 .seealso: `PetscDTIndexToBary` 3331 @*/ 3332 PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index) 3333 { 3334 PetscInt c; 3335 PetscInt i; 3336 PetscInt total; 3337 3338 PetscFunctionBeginHot; 3339 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 3340 if (!len) { 3341 PetscCheck(!sum, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 3342 *index = 0; 3343 PetscFunctionReturn(PETSC_SUCCESS); 3344 } 3345 for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c; 3346 i = total - 1; 3347 c = len - 1; 3348 sum -= coord[c]; 3349 while (sum > 0) { 3350 PetscInt subtotal; 3351 PetscInt s; 3352 3353 for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s; 3354 i -= subtotal; 3355 sum -= coord[--c]; 3356 } 3357 *index = i; 3358 PetscFunctionReturn(PETSC_SUCCESS); 3359 } 3360 3361 /*@ 3362 PetscQuadratureComputePermutations - Compute permutations of quadrature points corresponding to domain orientations 3363 3364 Input Parameter: 3365 . quad - The `PetscQuadrature` 3366 3367 Output Parameters: 3368 + Np - The number of domain orientations 3369 - perm - An array of `IS` permutations, one for ech orientation, 3370 3371 Level: developer 3372 3373 .seealso: `PetscQuadratureSetCellType()`, `PetscQuadrature` 3374 @*/ 3375 PetscErrorCode PetscQuadratureComputePermutations(PetscQuadrature quad, PeOp PetscInt *Np, IS *perm[]) 3376 { 3377 DMPolytopeType ct; 3378 const PetscReal *xq, *wq; 3379 PetscInt dim, qdim, d, Na, o, Nq, q, qp; 3380 3381 PetscFunctionBegin; 3382 PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &xq, &wq)); 3383 PetscCall(PetscQuadratureGetCellType(quad, &ct)); 3384 dim = DMPolytopeTypeGetDim(ct); 3385 Na = DMPolytopeTypeGetNumArrangements(ct); 3386 PetscCall(PetscMalloc1(Na, perm)); 3387 if (Np) *Np = Na; 3388 Na /= 2; 3389 for (o = -Na; o < Na; ++o) { 3390 DM refdm; 3391 PetscInt *idx; 3392 PetscReal xi0[3] = {-1., -1., -1.}, v0[3], J[9], detJ, txq[3]; 3393 PetscBool flg; 3394 3395 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm)); 3396 PetscCall(DMPlexOrientPoint(refdm, 0, o)); 3397 PetscCall(DMPlexComputeCellGeometryFEM(refdm, 0, NULL, v0, J, NULL, &detJ)); 3398 PetscCall(PetscMalloc1(Nq, &idx)); 3399 for (q = 0; q < Nq; ++q) { 3400 CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q * dim], txq); 3401 for (qp = 0; qp < Nq; ++qp) { 3402 PetscReal diff = 0.; 3403 3404 for (d = 0; d < dim; ++d) diff += PetscAbsReal(txq[d] - xq[qp * dim + d]); 3405 if (diff < PETSC_SMALL) break; 3406 } 3407 PetscCheck(qp < Nq, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Transformed quad point %" PetscInt_FMT " does not match another quad point", q); 3408 idx[q] = qp; 3409 } 3410 PetscCall(DMDestroy(&refdm)); 3411 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nq, idx, PETSC_OWN_POINTER, &(*perm)[o + Na])); 3412 PetscCall(ISGetInfo((*perm)[o + Na], IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg)); 3413 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ordering for orientation %" PetscInt_FMT " was not a permutation", o); 3414 PetscCall(ISSetPermutation((*perm)[o + Na])); 3415 } 3416 if (!Na) (*perm)[0] = NULL; 3417 PetscFunctionReturn(PETSC_SUCCESS); 3418 } 3419 3420 /*@ 3421 PetscDTCreateQuadratureByCell - Create default quadrature for a given cell 3422 3423 Not collective 3424 3425 Input Parameters: 3426 + ct - The integration domain 3427 . qorder - The desired quadrature order 3428 - qtype - The type of simplex quadrature, or PETSCDTSIMPLEXQUAD_DEFAULT 3429 3430 Output Parameters: 3431 + q - The cell quadrature 3432 - fq - The face quadrature 3433 3434 Level: developer 3435 3436 .seealso: `PetscDTCreateDefaultQuadrature()`, `PetscFECreateDefault()`, `PetscDTGaussTensorQuadrature()`, `PetscDTSimplexQuadrature()`, `PetscDTTensorQuadratureCreate()` 3437 @*/ 3438 PetscErrorCode PetscDTCreateQuadratureByCell(DMPolytopeType ct, PetscInt qorder, PetscDTSimplexQuadratureType qtype, PetscQuadrature *q, PetscQuadrature *fq) 3439 { 3440 const PetscInt quadPointsPerEdge = PetscMax(qorder + 1, 1); 3441 const PetscInt dim = DMPolytopeTypeGetDim(ct); 3442 3443 PetscFunctionBegin; 3444 switch (ct) { 3445 case DM_POLYTOPE_SEGMENT: 3446 case DM_POLYTOPE_POINT_PRISM_TENSOR: 3447 case DM_POLYTOPE_QUADRILATERAL: 3448 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3449 case DM_POLYTOPE_HEXAHEDRON: 3450 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3451 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q)); 3452 PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq)); 3453 break; 3454 case DM_POLYTOPE_TRIANGLE: 3455 case DM_POLYTOPE_TETRAHEDRON: 3456 PetscCall(PetscDTSimplexQuadrature(dim, 2 * qorder, qtype, q)); 3457 PetscCall(PetscDTSimplexQuadrature(dim - 1, 2 * qorder, qtype, fq)); 3458 break; 3459 case DM_POLYTOPE_TRI_PRISM: 3460 case DM_POLYTOPE_TRI_PRISM_TENSOR: { 3461 PetscQuadrature q1, q2; 3462 3463 // TODO: this should be able to use symmetric rules, but doing so causes tests to fail 3464 PetscCall(PetscDTSimplexQuadrature(2, 2 * qorder, PETSCDTSIMPLEXQUAD_CONIC, &q1)); 3465 PetscCall(PetscDTGaussTensorQuadrature(1, 1, quadPointsPerEdge, -1.0, 1.0, &q2)); 3466 PetscCall(PetscDTTensorQuadratureCreate(q1, q2, q)); 3467 PetscCall(PetscQuadratureDestroy(&q2)); 3468 *fq = q1; 3469 /* TODO Need separate quadratures for each face */ 3470 } break; 3471 default: 3472 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No quadrature for celltype %s", DMPolytopeTypes[PetscMin(ct, DM_POLYTOPE_UNKNOWN)]); 3473 } 3474 PetscFunctionReturn(PETSC_SUCCESS); 3475 } 3476 3477 /*@ 3478 PetscDTCreateDefaultQuadrature - Create default quadrature for a given cell 3479 3480 Not collective 3481 3482 Input Parameters: 3483 + ct - The integration domain 3484 - qorder - The desired quadrature order 3485 3486 Output Parameters: 3487 + q - The cell quadrature 3488 - fq - The face quadrature 3489 3490 Level: developer 3491 3492 .seealso: `PetscDTCreateQuadratureByCell()`, `PetscFECreateDefault()`, `PetscDTGaussTensorQuadrature()`, `PetscDTSimplexQuadrature()`, `PetscDTTensorQuadratureCreate()` 3493 @*/ 3494 PetscErrorCode PetscDTCreateDefaultQuadrature(DMPolytopeType ct, PetscInt qorder, PetscQuadrature *q, PetscQuadrature *fq) 3495 { 3496 PetscFunctionBegin; 3497 PetscCall(PetscDTCreateQuadratureByCell(ct, qorder, PETSCDTSIMPLEXQUAD_DEFAULT, q, fq)); 3498 PetscFunctionReturn(PETSC_SUCCESS); 3499 } 3500