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