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 - Return 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 Notes: 300 From Fortran you must 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), it is assumed that they represent multiple k-forms) [see `PetscDTAltVPullback()` for interpretation of formDegree] 478 479 Output Parameter: 480 . 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 been pulled-back by the pseudoinverse of `J` to the k-form weights in the image space. 481 482 Level: intermediate 483 484 Note: 485 The new quadrature rule will have a different number of components if spaces have different dimensions. 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. 486 487 .seealso: `PetscQuadrature`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()` 488 @*/ 489 PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq) 490 { 491 PetscInt dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c; 492 const PetscReal *points; 493 const PetscReal *weights; 494 PetscReal *imagePoints, *imageWeights; 495 PetscReal *Jinv; 496 PetscReal *Jinvstar; 497 498 PetscFunctionBegin; 499 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 500 PetscCheck(imageDim >= PetscAbsInt(formDegree), PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Cannot represent a %" PetscInt_FMT "-form in %" PetscInt_FMT " dimensions", PetscAbsInt(formDegree), imageDim); 501 PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights)); 502 PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize)); 503 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); 504 Ncopies = Nc / formSize; 505 PetscCall(PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize)); 506 imageNc = Ncopies * imageFormSize; 507 PetscCall(PetscMalloc1(Npoints * imageDim, &imagePoints)); 508 PetscCall(PetscMalloc1(Npoints * imageNc, &imageWeights)); 509 PetscCall(PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar)); 510 PetscCall(PetscDTJacobianInverse_Internal(imageDim, dim, J, Jinv)); 511 PetscCall(PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar)); 512 for (pt = 0; pt < Npoints; pt++) { 513 const PetscReal *point = &points[pt * dim]; 514 PetscReal *imagePoint = &imagePoints[pt * imageDim]; 515 516 for (i = 0; i < imageDim; i++) { 517 PetscReal val = originImage[i]; 518 519 for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]); 520 imagePoint[i] = val; 521 } 522 for (c = 0; c < Ncopies; c++) { 523 const PetscReal *form = &weights[pt * Nc + c * formSize]; 524 PetscReal *imageForm = &imageWeights[pt * imageNc + c * imageFormSize]; 525 526 for (i = 0; i < imageFormSize; i++) { 527 PetscReal val = 0.; 528 529 for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j]; 530 imageForm[i] = val; 531 } 532 } 533 } 534 PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq)); 535 PetscCall(PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights)); 536 PetscCall(PetscFree2(Jinv, Jinvstar)); 537 PetscFunctionReturn(PETSC_SUCCESS); 538 } 539 540 /*@C 541 PetscQuadratureSetData - Sets the data defining the quadrature 542 543 Not Collective 544 545 Input Parameters: 546 + q - The `PetscQuadrature` object 547 . dim - The spatial dimension 548 . Nc - The number of components 549 . npoints - The number of quadrature points 550 . points - The coordinates of each quadrature point 551 - weights - The weight of each quadrature point 552 553 Level: intermediate 554 555 Note: 556 This routine owns the references to points and weights, so they must be allocated using `PetscMalloc()` and the user should not free them. 557 558 .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()` 559 @*/ 560 PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[]) 561 { 562 PetscFunctionBegin; 563 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 564 if (dim >= 0) q->dim = dim; 565 if (Nc >= 0) q->Nc = Nc; 566 if (npoints >= 0) q->numPoints = npoints; 567 if (points) { 568 PetscAssertPointer(points, 5); 569 q->points = points; 570 } 571 if (weights) { 572 PetscAssertPointer(weights, 6); 573 q->weights = weights; 574 } 575 PetscFunctionReturn(PETSC_SUCCESS); 576 } 577 578 static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v) 579 { 580 PetscInt q, d, c; 581 PetscViewerFormat format; 582 583 PetscFunctionBegin; 584 if (quad->Nc > 1) 585 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)); 586 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)); 587 PetscCall(PetscViewerGetFormat(v, &format)); 588 if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS); 589 for (q = 0; q < quad->numPoints; ++q) { 590 PetscCall(PetscViewerASCIIPrintf(v, "p%" PetscInt_FMT " (", q)); 591 PetscCall(PetscViewerASCIIUseTabs(v, PETSC_FALSE)); 592 for (d = 0; d < quad->dim; ++d) { 593 if (d) PetscCall(PetscViewerASCIIPrintf(v, ", ")); 594 PetscCall(PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q * quad->dim + d])); 595 } 596 PetscCall(PetscViewerASCIIPrintf(v, ") ")); 597 if (quad->Nc > 1) PetscCall(PetscViewerASCIIPrintf(v, "w%" PetscInt_FMT " (", q)); 598 for (c = 0; c < quad->Nc; ++c) { 599 if (c) PetscCall(PetscViewerASCIIPrintf(v, ", ")); 600 PetscCall(PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q * quad->Nc + c])); 601 } 602 if (quad->Nc > 1) PetscCall(PetscViewerASCIIPrintf(v, ")")); 603 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 604 PetscCall(PetscViewerASCIIUseTabs(v, PETSC_TRUE)); 605 } 606 PetscFunctionReturn(PETSC_SUCCESS); 607 } 608 609 /*@C 610 PetscQuadratureView - View a `PetscQuadrature` object 611 612 Collective 613 614 Input Parameters: 615 + quad - The `PetscQuadrature` object 616 - viewer - The `PetscViewer` object 617 618 Level: beginner 619 620 .seealso: `PetscQuadrature`, `PetscViewer`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()` 621 @*/ 622 PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer) 623 { 624 PetscBool iascii; 625 626 PetscFunctionBegin; 627 PetscValidHeader(quad, 1); 628 if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 629 if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)quad), &viewer)); 630 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 631 PetscCall(PetscViewerASCIIPushTab(viewer)); 632 if (iascii) PetscCall(PetscQuadratureView_Ascii(quad, viewer)); 633 PetscCall(PetscViewerASCIIPopTab(viewer)); 634 PetscFunctionReturn(PETSC_SUCCESS); 635 } 636 637 /*@C 638 PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement 639 640 Not Collective; No Fortran Support 641 642 Input Parameters: 643 + q - The original `PetscQuadrature` 644 . numSubelements - The number of subelements the original element is divided into 645 . v0 - An array of the initial points for each subelement 646 - jac - An array of the Jacobian mappings from the reference to each subelement 647 648 Output Parameter: 649 . qref - The dimension 650 651 Level: intermediate 652 653 Note: 654 Together v0 and jac define an affine mapping from the original reference element to each subelement 655 656 .seealso: `PetscQuadrature`, `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()` 657 @*/ 658 PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref) 659 { 660 DMPolytopeType ct; 661 const PetscReal *points, *weights; 662 PetscReal *pointsRef, *weightsRef; 663 PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e; 664 665 PetscFunctionBegin; 666 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 667 PetscAssertPointer(v0, 3); 668 PetscAssertPointer(jac, 4); 669 PetscAssertPointer(qref, 5); 670 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, qref)); 671 PetscCall(PetscQuadratureGetCellType(q, &ct)); 672 PetscCall(PetscQuadratureGetOrder(q, &order)); 673 PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights)); 674 npointsRef = npoints * numSubelements; 675 PetscCall(PetscMalloc1(npointsRef * dim, &pointsRef)); 676 PetscCall(PetscMalloc1(npointsRef * Nc, &weightsRef)); 677 for (c = 0; c < numSubelements; ++c) { 678 for (p = 0; p < npoints; ++p) { 679 for (d = 0; d < dim; ++d) { 680 pointsRef[(c * npoints + p) * dim + d] = v0[c * dim + d]; 681 for (e = 0; e < dim; ++e) pointsRef[(c * npoints + p) * dim + d] += jac[(c * dim + d) * dim + e] * (points[p * dim + e] + 1.0); 682 } 683 /* Could also use detJ here */ 684 for (cp = 0; cp < Nc; ++cp) weightsRef[(c * npoints + p) * Nc + cp] = weights[p * Nc + cp] / numSubelements; 685 } 686 } 687 PetscCall(PetscQuadratureSetCellType(*qref, ct)); 688 PetscCall(PetscQuadratureSetOrder(*qref, order)); 689 PetscCall(PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef)); 690 PetscFunctionReturn(PETSC_SUCCESS); 691 } 692 693 /* Compute the coefficients for the Jacobi polynomial recurrence, 694 * 695 * J^{a,b}_n(x) = (cnm1 + cnm1x * x) * J^{a,b}_{n-1}(x) - cnm2 * J^{a,b}_{n-2}(x). 696 */ 697 #define PetscDTJacobiRecurrence_Internal(n, a, b, cnm1, cnm1x, cnm2) \ 698 do { \ 699 PetscReal _a = (a); \ 700 PetscReal _b = (b); \ 701 PetscReal _n = (n); \ 702 if (n == 1) { \ 703 (cnm1) = (_a - _b) * 0.5; \ 704 (cnm1x) = (_a + _b + 2.) * 0.5; \ 705 (cnm2) = 0.; \ 706 } else { \ 707 PetscReal _2n = _n + _n; \ 708 PetscReal _d = (_2n * (_n + _a + _b) * (_2n + _a + _b - 2)); \ 709 PetscReal _n1 = (_2n + _a + _b - 1.) * (_a * _a - _b * _b); \ 710 PetscReal _n1x = (_2n + _a + _b - 1.) * (_2n + _a + _b) * (_2n + _a + _b - 2); \ 711 PetscReal _n2 = 2. * ((_n + _a - 1.) * (_n + _b - 1.) * (_2n + _a + _b)); \ 712 (cnm1) = _n1 / _d; \ 713 (cnm1x) = _n1x / _d; \ 714 (cnm2) = _n2 / _d; \ 715 } \ 716 } while (0) 717 718 /*@ 719 PetscDTJacobiNorm - Compute the weighted L2 norm of a Jacobi polynomial. 720 721 $\| P^{\alpha,\beta}_n \|_{\alpha,\beta}^2 = \int_{-1}^1 (1 + x)^{\alpha} (1 - x)^{\beta} P^{\alpha,\beta}_n (x)^2 dx.$ 722 723 Input Parameters: 724 + alpha - the left exponent > -1 725 . beta - the right exponent > -1 726 - n - the polynomial degree 727 728 Output Parameter: 729 . norm - the weighted L2 norm 730 731 Level: beginner 732 733 .seealso: `PetscQuadrature`, `PetscDTJacobiEval()` 734 @*/ 735 PetscErrorCode PetscDTJacobiNorm(PetscReal alpha, PetscReal beta, PetscInt n, PetscReal *norm) 736 { 737 PetscReal twoab1; 738 PetscReal gr; 739 740 PetscFunctionBegin; 741 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent alpha %g <= -1. invalid", (double)alpha); 742 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent beta %g <= -1. invalid", (double)beta); 743 PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "n %" PetscInt_FMT " < 0 invalid", n); 744 twoab1 = PetscPowReal(2., alpha + beta + 1.); 745 #if defined(PETSC_HAVE_LGAMMA) 746 if (!n) { 747 gr = PetscExpReal(PetscLGamma(alpha + 1.) + PetscLGamma(beta + 1.) - PetscLGamma(alpha + beta + 2.)); 748 } else { 749 gr = PetscExpReal(PetscLGamma(n + alpha + 1.) + PetscLGamma(n + beta + 1.) - (PetscLGamma(n + 1.) + PetscLGamma(n + alpha + beta + 1.))) / (n + n + alpha + beta + 1.); 750 } 751 #else 752 { 753 PetscInt alphai = (PetscInt)alpha; 754 PetscInt betai = (PetscInt)beta; 755 PetscInt i; 756 757 gr = n ? (1. / (n + n + alpha + beta + 1.)) : 1.; 758 if ((PetscReal)alphai == alpha) { 759 if (!n) { 760 for (i = 0; i < alphai; i++) gr *= (i + 1.) / (beta + i + 1.); 761 gr /= (alpha + beta + 1.); 762 } else { 763 for (i = 0; i < alphai; i++) gr *= (n + i + 1.) / (n + beta + i + 1.); 764 } 765 } else if ((PetscReal)betai == beta) { 766 if (!n) { 767 for (i = 0; i < betai; i++) gr *= (i + 1.) / (alpha + i + 2.); 768 gr /= (alpha + beta + 1.); 769 } else { 770 for (i = 0; i < betai; i++) gr *= (n + i + 1.) / (n + alpha + i + 1.); 771 } 772 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable."); 773 } 774 #endif 775 *norm = PetscSqrtReal(twoab1 * gr); 776 PetscFunctionReturn(PETSC_SUCCESS); 777 } 778 779 static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p) 780 { 781 PetscReal ak, bk; 782 PetscReal abk1; 783 PetscInt i, l, maxdegree; 784 785 PetscFunctionBegin; 786 maxdegree = degrees[ndegree - 1] - k; 787 ak = a + k; 788 bk = b + k; 789 abk1 = a + b + k + 1.; 790 if (maxdegree < 0) { 791 for (i = 0; i < npoints; i++) 792 for (l = 0; l < ndegree; l++) p[i * ndegree + l] = 0.; 793 PetscFunctionReturn(PETSC_SUCCESS); 794 } 795 for (i = 0; i < npoints; i++) { 796 PetscReal pm1, pm2, x; 797 PetscReal cnm1, cnm1x, cnm2; 798 PetscInt j, m; 799 800 x = points[i]; 801 pm2 = 1.; 802 PetscDTJacobiRecurrence_Internal(1, ak, bk, cnm1, cnm1x, cnm2); 803 pm1 = (cnm1 + cnm1x * x); 804 l = 0; 805 while (l < ndegree && degrees[l] - k < 0) p[l++] = 0.; 806 while (l < ndegree && degrees[l] - k == 0) { 807 p[l] = pm2; 808 for (m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5; 809 l++; 810 } 811 while (l < ndegree && degrees[l] - k == 1) { 812 p[l] = pm1; 813 for (m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5; 814 l++; 815 } 816 for (j = 2; j <= maxdegree; j++) { 817 PetscReal pp; 818 819 PetscDTJacobiRecurrence_Internal(j, ak, bk, cnm1, cnm1x, cnm2); 820 pp = (cnm1 + cnm1x * x) * pm1 - cnm2 * pm2; 821 pm2 = pm1; 822 pm1 = pp; 823 while (l < ndegree && degrees[l] - k == j) { 824 p[l] = pp; 825 for (m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5; 826 l++; 827 } 828 } 829 p += ndegree; 830 } 831 PetscFunctionReturn(PETSC_SUCCESS); 832 } 833 834 /*@ 835 PetscDTJacobiEvalJet - Evaluate the jet (function and derivatives) of the Jacobi polynomials basis up to a given degree. 836 837 Input Parameters: 838 + alpha - the left exponent of the weight 839 . beta - the right exponetn of the weight 840 . npoints - the number of points to evaluate the polynomials at 841 . points - [npoints] array of point coordinates 842 . degree - the maximm degree polynomial space to evaluate, (degree + 1) will be evaluated total. 843 - k - the maximum derivative to evaluate in the jet, (k + 1) will be evaluated total. 844 845 Output Parameter: 846 . p - an array containing the evaluations of the Jacobi polynomials's jets on the points. the size is (degree + 1) x 847 (k + 1) x npoints, which also describes the order of the dimensions of this three-dimensional array: the first 848 (slowest varying) dimension is polynomial degree; the second dimension is derivative order; the third (fastest 849 varying) dimension is the index of the evaluation point. 850 851 Level: advanced 852 853 Notes: 854 The Jacobi polynomials with indices $\alpha$ and $\beta$ are orthogonal with respect to the 855 weighted inner product $\langle f, g \rangle = \int_{-1}^1 (1+x)^{\alpha} (1-x)^{\beta} f(x) 856 g(x) dx$. 857 858 .seealso: `PetscDTJacobiEval()`, `PetscDTPKDEvalJet()` 859 @*/ 860 PetscErrorCode PetscDTJacobiEvalJet(PetscReal alpha, PetscReal beta, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[]) 861 { 862 PetscInt i, j, l; 863 PetscInt *degrees; 864 PetscReal *psingle; 865 866 PetscFunctionBegin; 867 if (degree == 0) { 868 PetscInt zero = 0; 869 870 for (i = 0; i <= k; i++) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, 1, &zero, &p[i * npoints])); 871 PetscFunctionReturn(PETSC_SUCCESS); 872 } 873 PetscCall(PetscMalloc1(degree + 1, °rees)); 874 PetscCall(PetscMalloc1((degree + 1) * npoints, &psingle)); 875 for (i = 0; i <= degree; i++) degrees[i] = i; 876 for (i = 0; i <= k; i++) { 877 PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, degree + 1, degrees, psingle)); 878 for (j = 0; j <= degree; j++) { 879 for (l = 0; l < npoints; l++) p[(j * (k + 1) + i) * npoints + l] = psingle[l * (degree + 1) + j]; 880 } 881 } 882 PetscCall(PetscFree(psingle)); 883 PetscCall(PetscFree(degrees)); 884 PetscFunctionReturn(PETSC_SUCCESS); 885 } 886 887 /*@ 888 PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$ at a set of points 889 at points 890 891 Not Collective 892 893 Input Parameters: 894 + npoints - number of spatial points to evaluate at 895 . alpha - the left exponent > -1 896 . beta - the right exponent > -1 897 . points - array of locations to evaluate at 898 . ndegree - number of basis degrees to evaluate 899 - degrees - sorted array of degrees to evaluate 900 901 Output Parameters: 902 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL) 903 . D - row-oriented derivative evaluation matrix (or NULL) 904 - D2 - row-oriented second derivative evaluation matrix (or NULL) 905 906 Level: intermediate 907 908 .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()` 909 @*/ 910 PetscErrorCode PetscDTJacobiEval(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *B, PetscReal *D, PetscReal *D2) 911 { 912 PetscFunctionBegin; 913 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1."); 914 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1."); 915 if (!npoints || !ndegree) PetscFunctionReturn(PETSC_SUCCESS); 916 if (B) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B)); 917 if (D) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D)); 918 if (D2) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2)); 919 PetscFunctionReturn(PETSC_SUCCESS); 920 } 921 922 /*@ 923 PetscDTLegendreEval - evaluate Legendre polynomials at points 924 925 Not Collective 926 927 Input Parameters: 928 + npoints - number of spatial points to evaluate at 929 . points - array of locations to evaluate at 930 . ndegree - number of basis degrees to evaluate 931 - degrees - sorted array of degrees to evaluate 932 933 Output Parameters: 934 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL) 935 . D - row-oriented derivative evaluation matrix (or NULL) 936 - D2 - row-oriented second derivative evaluation matrix (or NULL) 937 938 Level: intermediate 939 940 .seealso: `PetscDTGaussQuadrature()` 941 @*/ 942 PetscErrorCode PetscDTLegendreEval(PetscInt npoints, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *B, PetscReal *D, PetscReal *D2) 943 { 944 PetscFunctionBegin; 945 PetscCall(PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2)); 946 PetscFunctionReturn(PETSC_SUCCESS); 947 } 948 949 /*@ 950 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, then the index of x is smaller than the index of y) 951 952 Input Parameters: 953 + len - the desired length of the degree tuple 954 - index - the index to convert: should be >= 0 955 956 Output Parameter: 957 . degtup - will be filled with a tuple of degrees 958 959 Level: beginner 960 961 Note: 962 For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples 963 acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the 964 last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1). 965 966 .seealso: `PetscDTGradedOrderToIndex()` 967 @*/ 968 PetscErrorCode PetscDTIndexToGradedOrder(PetscInt len, PetscInt index, PetscInt degtup[]) 969 { 970 PetscInt i, total; 971 PetscInt sum; 972 973 PetscFunctionBeginHot; 974 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 975 PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative"); 976 total = 1; 977 sum = 0; 978 while (index >= total) { 979 index -= total; 980 total = (total * (len + sum)) / (sum + 1); 981 sum++; 982 } 983 for (i = 0; i < len; i++) { 984 PetscInt c; 985 986 degtup[i] = sum; 987 for (c = 0, total = 1; c < sum; c++) { 988 /* going into the loop, total is the number of way to have a tuple of sum exactly c with length len - 1 - i */ 989 if (index < total) break; 990 index -= total; 991 total = (total * (len - 1 - i + c)) / (c + 1); 992 degtup[i]--; 993 } 994 sum -= degtup[i]; 995 } 996 PetscFunctionReturn(PETSC_SUCCESS); 997 } 998 999 /*@ 1000 PetscDTGradedOrderToIndex - convert a tuple into an index in a graded order, the inverse of `PetscDTIndexToGradedOrder()`. 1001 1002 Input Parameters: 1003 + len - the length of the degree tuple 1004 - degtup - tuple with this length 1005 1006 Output Parameter: 1007 . index - index in graded order: >= 0 1008 1009 Level: beginner 1010 1011 Note: 1012 For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples 1013 acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the 1014 last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1). 1015 1016 .seealso: `PetscDTIndexToGradedOrder()` 1017 @*/ 1018 PetscErrorCode PetscDTGradedOrderToIndex(PetscInt len, const PetscInt degtup[], PetscInt *index) 1019 { 1020 PetscInt i, idx, sum, total; 1021 1022 PetscFunctionBeginHot; 1023 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 1024 for (i = 0, sum = 0; i < len; i++) sum += degtup[i]; 1025 idx = 0; 1026 total = 1; 1027 for (i = 0; i < sum; i++) { 1028 idx += total; 1029 total = (total * (len + i)) / (i + 1); 1030 } 1031 for (i = 0; i < len - 1; i++) { 1032 PetscInt c; 1033 1034 total = 1; 1035 sum -= degtup[i]; 1036 for (c = 0; c < sum; c++) { 1037 idx += total; 1038 total = (total * (len - 1 - i + c)) / (c + 1); 1039 } 1040 } 1041 *index = idx; 1042 PetscFunctionReturn(PETSC_SUCCESS); 1043 } 1044 1045 static PetscBool PKDCite = PETSC_FALSE; 1046 const char PKDCitation[] = "@article{Kirby2010,\n" 1047 " title={Singularity-free evaluation of collapsed-coordinate orthogonal polynomials},\n" 1048 " author={Kirby, Robert C},\n" 1049 " journal={ACM Transactions on Mathematical Software (TOMS)},\n" 1050 " volume={37},\n" 1051 " number={1},\n" 1052 " pages={1--16},\n" 1053 " year={2010},\n" 1054 " publisher={ACM New York, NY, USA}\n}\n"; 1055 1056 /*@ 1057 PetscDTPKDEvalJet - Evaluate the jet (function and derivatives) of the Proriol-Koornwinder-Dubiner (PKD) basis for 1058 the space of polynomials up to a given degree. 1059 1060 Input Parameters: 1061 + dim - the number of variables in the multivariate polynomials 1062 . npoints - the number of points to evaluate the polynomials at 1063 . points - [npoints x dim] array of point coordinates 1064 . 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. 1065 - k - the maximum order partial derivative to evaluate in the jet. There are (dim + k choose dim) partial derivatives 1066 in the jet. Choosing k = 0 means to evaluate just the function and no derivatives 1067 1068 Output Parameter: 1069 . p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is ((dim + degree) 1070 choose dim) x ((dim + k) choose dim) x npoints, which also describes the order of the dimensions of this 1071 three-dimensional array: the first (slowest varying) dimension is basis function index; the second dimension is jet 1072 index; the third (fastest varying) dimension is the index of the evaluation point. 1073 1074 Level: advanced 1075 1076 Notes: 1077 The PKD basis is L2-orthonormal on the biunit simplex (which is used as the reference element 1078 for finite elements in PETSc), which makes it a stable basis to use for evaluating 1079 polynomials in that domain. 1080 1081 The ordering of the basis functions, and the ordering of the derivatives in the jet, both follow the graded 1082 ordering of `PetscDTIndexToGradedOrder()` and `PetscDTGradedOrderToIndex()`. For example, in 3D, the polynomial with 1083 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); 1084 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). 1085 1086 The implementation uses Kirby's singularity-free evaluation algorithm, https://doi.org/10.1145/1644001.1644006. 1087 1088 .seealso: `PetscDTGradedOrderToIndex()`, `PetscDTIndexToGradedOrder()`, `PetscDTJacobiEvalJet()` 1089 @*/ 1090 PetscErrorCode PetscDTPKDEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[]) 1091 { 1092 PetscInt degidx, kidx, d, pt; 1093 PetscInt Nk, Ndeg; 1094 PetscInt *ktup, *degtup; 1095 PetscReal *scales, initscale, scaleexp; 1096 1097 PetscFunctionBegin; 1098 PetscCall(PetscCitationsRegister(PKDCitation, &PKDCite)); 1099 PetscCall(PetscDTBinomialInt(dim + k, k, &Nk)); 1100 PetscCall(PetscDTBinomialInt(degree + dim, degree, &Ndeg)); 1101 PetscCall(PetscMalloc2(dim, °tup, dim, &ktup)); 1102 PetscCall(PetscMalloc1(Ndeg, &scales)); 1103 initscale = 1.; 1104 if (dim > 1) { 1105 PetscCall(PetscDTBinomial(dim, 2, &scaleexp)); 1106 initscale = PetscPowReal(2., scaleexp * 0.5); 1107 } 1108 for (degidx = 0; degidx < Ndeg; degidx++) { 1109 PetscInt e, i; 1110 PetscInt m1idx = -1, m2idx = -1; 1111 PetscInt n; 1112 PetscInt degsum; 1113 PetscReal alpha; 1114 PetscReal cnm1, cnm1x, cnm2; 1115 PetscReal norm; 1116 1117 PetscCall(PetscDTIndexToGradedOrder(dim, degidx, degtup)); 1118 for (d = dim - 1; d >= 0; d--) 1119 if (degtup[d]) break; 1120 if (d < 0) { /* constant is 1 everywhere, all derivatives are zero */ 1121 scales[degidx] = initscale; 1122 for (e = 0; e < dim; e++) { 1123 PetscCall(PetscDTJacobiNorm(e, 0., 0, &norm)); 1124 scales[degidx] /= norm; 1125 } 1126 for (i = 0; i < npoints; i++) p[degidx * Nk * npoints + i] = 1.; 1127 for (i = 0; i < (Nk - 1) * npoints; i++) p[(degidx * Nk + 1) * npoints + i] = 0.; 1128 continue; 1129 } 1130 n = degtup[d]; 1131 degtup[d]--; 1132 PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m1idx)); 1133 if (degtup[d] > 0) { 1134 degtup[d]--; 1135 PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m2idx)); 1136 degtup[d]++; 1137 } 1138 degtup[d]++; 1139 for (e = 0, degsum = 0; e < d; e++) degsum += degtup[e]; 1140 alpha = 2 * degsum + d; 1141 PetscDTJacobiRecurrence_Internal(n, alpha, 0., cnm1, cnm1x, cnm2); 1142 1143 scales[degidx] = initscale; 1144 for (e = 0, degsum = 0; e < dim; e++) { 1145 PetscInt f; 1146 PetscReal ealpha; 1147 PetscReal enorm; 1148 1149 ealpha = 2 * degsum + e; 1150 for (f = 0; f < degsum; f++) scales[degidx] *= 2.; 1151 PetscCall(PetscDTJacobiNorm(ealpha, 0., degtup[e], &enorm)); 1152 scales[degidx] /= enorm; 1153 degsum += degtup[e]; 1154 } 1155 1156 for (pt = 0; pt < npoints; pt++) { 1157 /* compute the multipliers */ 1158 PetscReal thetanm1, thetanm1x, thetanm2; 1159 1160 thetanm1x = dim - (d + 1) + 2. * points[pt * dim + d]; 1161 for (e = d + 1; e < dim; e++) thetanm1x += points[pt * dim + e]; 1162 thetanm1x *= 0.5; 1163 thetanm1 = (2. - (dim - (d + 1))); 1164 for (e = d + 1; e < dim; e++) thetanm1 -= points[pt * dim + e]; 1165 thetanm1 *= 0.5; 1166 thetanm2 = thetanm1 * thetanm1; 1167 1168 for (kidx = 0; kidx < Nk; kidx++) { 1169 PetscInt f; 1170 1171 PetscCall(PetscDTIndexToGradedOrder(dim, kidx, ktup)); 1172 /* first sum in the same derivative terms */ 1173 p[(degidx * Nk + kidx) * npoints + pt] = (cnm1 * thetanm1 + cnm1x * thetanm1x) * p[(m1idx * Nk + kidx) * npoints + pt]; 1174 if (m2idx >= 0) p[(degidx * Nk + kidx) * npoints + pt] -= cnm2 * thetanm2 * p[(m2idx * Nk + kidx) * npoints + pt]; 1175 1176 for (f = d; f < dim; f++) { 1177 PetscInt km1idx, mplty = ktup[f]; 1178 1179 if (!mplty) continue; 1180 ktup[f]--; 1181 PetscCall(PetscDTGradedOrderToIndex(dim, ktup, &km1idx)); 1182 1183 /* the derivative of cnm1x * thetanm1x wrt x variable f is 0.5 * cnm1x if f > d otherwise it is cnm1x */ 1184 /* the derivative of cnm1 * thetanm1 wrt x variable f is 0 if f == d, otherwise it is -0.5 * cnm1 */ 1185 /* the derivative of -cnm2 * thetanm2 wrt x variable f is 0 if f == d, otherwise it is cnm2 * thetanm1 */ 1186 if (f > d) { 1187 PetscInt f2; 1188 1189 p[(degidx * Nk + kidx) * npoints + pt] += mplty * 0.5 * (cnm1x - cnm1) * p[(m1idx * Nk + km1idx) * npoints + pt]; 1190 if (m2idx >= 0) { 1191 p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm2 * thetanm1 * p[(m2idx * Nk + km1idx) * npoints + pt]; 1192 /* second derivatives of -cnm2 * thetanm2 wrt x variable f,f2 is like - 0.5 * cnm2 */ 1193 for (f2 = f; f2 < dim; f2++) { 1194 PetscInt km2idx, mplty2 = ktup[f2]; 1195 PetscInt factor; 1196 1197 if (!mplty2) continue; 1198 ktup[f2]--; 1199 PetscCall(PetscDTGradedOrderToIndex(dim, ktup, &km2idx)); 1200 1201 factor = mplty * mplty2; 1202 if (f == f2) factor /= 2; 1203 p[(degidx * Nk + kidx) * npoints + pt] -= 0.5 * factor * cnm2 * p[(m2idx * Nk + km2idx) * npoints + pt]; 1204 ktup[f2]++; 1205 } 1206 } 1207 } else { 1208 p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm1x * p[(m1idx * Nk + km1idx) * npoints + pt]; 1209 } 1210 ktup[f]++; 1211 } 1212 } 1213 } 1214 } 1215 for (degidx = 0; degidx < Ndeg; degidx++) { 1216 PetscReal scale = scales[degidx]; 1217 PetscInt i; 1218 1219 for (i = 0; i < Nk * npoints; i++) p[degidx * Nk * npoints + i] *= scale; 1220 } 1221 PetscCall(PetscFree(scales)); 1222 PetscCall(PetscFree2(degtup, ktup)); 1223 PetscFunctionReturn(PETSC_SUCCESS); 1224 } 1225 1226 /*@ 1227 PetscDTPTrimmedSize - The size of the trimmed polynomial space of k-forms with a given degree and form degree, 1228 which can be evaluated in `PetscDTPTrimmedEvalJet()`. 1229 1230 Input Parameters: 1231 + dim - the number of variables in the multivariate polynomials 1232 . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space. 1233 - formDegree - the degree of the form 1234 1235 Output Parameter: 1236 . size - The number ((`dim` + `degree`) choose (`dim` + `formDegree`)) x ((`degree` + `formDegree` - 1) choose (`formDegree`)) 1237 1238 Level: advanced 1239 1240 .seealso: `PetscDTPTrimmedEvalJet()` 1241 @*/ 1242 PetscErrorCode PetscDTPTrimmedSize(PetscInt dim, PetscInt degree, PetscInt formDegree, PetscInt *size) 1243 { 1244 PetscInt Nrk, Nbpt; // number of trimmed polynomials 1245 1246 PetscFunctionBegin; 1247 formDegree = PetscAbsInt(formDegree); 1248 PetscCall(PetscDTBinomialInt(degree + dim, degree + formDegree, &Nbpt)); 1249 PetscCall(PetscDTBinomialInt(degree + formDegree - 1, formDegree, &Nrk)); 1250 Nbpt *= Nrk; 1251 *size = Nbpt; 1252 PetscFunctionReturn(PETSC_SUCCESS); 1253 } 1254 1255 /* there was a reference implementation based on section 4.4 of Arnold, Falk & Winther (acta numerica, 2006), but it 1256 * was inferior to this implementation */ 1257 static PetscErrorCode PetscDTPTrimmedEvalJet_Internal(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[]) 1258 { 1259 PetscInt formDegreeOrig = formDegree; 1260 PetscBool formNegative = (formDegreeOrig < 0) ? PETSC_TRUE : PETSC_FALSE; 1261 1262 PetscFunctionBegin; 1263 formDegree = PetscAbsInt(formDegreeOrig); 1264 if (formDegree == 0) { 1265 PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree, jetDegree, p)); 1266 PetscFunctionReturn(PETSC_SUCCESS); 1267 } 1268 if (formDegree == dim) { 1269 PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p)); 1270 PetscFunctionReturn(PETSC_SUCCESS); 1271 } 1272 PetscInt Nbpt; 1273 PetscCall(PetscDTPTrimmedSize(dim, degree, formDegree, &Nbpt)); 1274 PetscInt Nf; 1275 PetscCall(PetscDTBinomialInt(dim, formDegree, &Nf)); 1276 PetscInt Nk; 1277 PetscCall(PetscDTBinomialInt(dim + jetDegree, dim, &Nk)); 1278 PetscCall(PetscArrayzero(p, Nbpt * Nf * Nk * npoints)); 1279 1280 PetscInt Nbpm1; // number of scalar polynomials up to degree - 1; 1281 PetscCall(PetscDTBinomialInt(dim + degree - 1, dim, &Nbpm1)); 1282 PetscReal *p_scalar; 1283 PetscCall(PetscMalloc1(Nbpm1 * Nk * npoints, &p_scalar)); 1284 PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p_scalar)); 1285 PetscInt total = 0; 1286 // First add the full polynomials up to degree - 1 into the basis: take the scalar 1287 // and copy one for each form component 1288 for (PetscInt i = 0; i < Nbpm1; i++) { 1289 const PetscReal *src = &p_scalar[i * Nk * npoints]; 1290 for (PetscInt f = 0; f < Nf; f++) { 1291 PetscReal *dest = &p[(total++ * Nf + f) * Nk * npoints]; 1292 PetscCall(PetscArraycpy(dest, src, Nk * npoints)); 1293 } 1294 } 1295 PetscInt *form_atoms; 1296 PetscCall(PetscMalloc1(formDegree + 1, &form_atoms)); 1297 // construct the interior product pattern 1298 PetscInt(*pattern)[3]; 1299 PetscInt Nf1; // number of formDegree + 1 forms 1300 PetscCall(PetscDTBinomialInt(dim, formDegree + 1, &Nf1)); 1301 PetscInt nnz = Nf1 * (formDegree + 1); 1302 PetscCall(PetscMalloc1(Nf1 * (formDegree + 1), &pattern)); 1303 PetscCall(PetscDTAltVInteriorPattern(dim, formDegree + 1, pattern)); 1304 PetscReal centroid = (1. - dim) / (dim + 1.); 1305 PetscInt *deriv; 1306 PetscCall(PetscMalloc1(dim, &deriv)); 1307 for (PetscInt d = dim; d >= formDegree + 1; d--) { 1308 PetscInt Nfd1; // number of formDegree + 1 forms in dimension d that include dx_0 1309 // (equal to the number of formDegree forms in dimension d-1) 1310 PetscCall(PetscDTBinomialInt(d - 1, formDegree, &Nfd1)); 1311 // The number of homogeneous (degree-1) scalar polynomials in d variables 1312 PetscInt Nh; 1313 PetscCall(PetscDTBinomialInt(d - 1 + degree - 1, d - 1, &Nh)); 1314 const PetscReal *h_scalar = &p_scalar[(Nbpm1 - Nh) * Nk * npoints]; 1315 for (PetscInt b = 0; b < Nh; b++) { 1316 const PetscReal *h_s = &h_scalar[b * Nk * npoints]; 1317 for (PetscInt f = 0; f < Nfd1; f++) { 1318 // construct all formDegree+1 forms that start with dx_(dim - d) /\ ... 1319 form_atoms[0] = dim - d; 1320 PetscCall(PetscDTEnumSubset(d - 1, formDegree, f, &form_atoms[1])); 1321 for (PetscInt i = 0; i < formDegree; i++) form_atoms[1 + i] += form_atoms[0] + 1; 1322 PetscInt f_ind; // index of the resulting form 1323 PetscCall(PetscDTSubsetIndex(dim, formDegree + 1, form_atoms, &f_ind)); 1324 PetscReal *p_f = &p[total++ * Nf * Nk * npoints]; 1325 for (PetscInt nz = 0; nz < nnz; nz++) { 1326 PetscInt i = pattern[nz][0]; // formDegree component 1327 PetscInt j = pattern[nz][1]; // (formDegree + 1) component 1328 PetscInt v = pattern[nz][2]; // coordinate component 1329 PetscReal scale = v < 0 ? -1. : 1.; 1330 1331 i = formNegative ? (Nf - 1 - i) : i; 1332 scale = (formNegative && (i & 1)) ? -scale : scale; 1333 v = v < 0 ? -(v + 1) : v; 1334 if (j != f_ind) continue; 1335 PetscReal *p_i = &p_f[i * Nk * npoints]; 1336 for (PetscInt jet = 0; jet < Nk; jet++) { 1337 const PetscReal *h_jet = &h_s[jet * npoints]; 1338 PetscReal *p_jet = &p_i[jet * npoints]; 1339 1340 for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * h_jet[pt] * (points[pt * dim + v] - centroid); 1341 PetscCall(PetscDTIndexToGradedOrder(dim, jet, deriv)); 1342 deriv[v]++; 1343 PetscReal mult = deriv[v]; 1344 PetscInt l; 1345 PetscCall(PetscDTGradedOrderToIndex(dim, deriv, &l)); 1346 if (l >= Nk) continue; 1347 p_jet = &p_i[l * npoints]; 1348 for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * mult * h_jet[pt]; 1349 deriv[v]--; 1350 } 1351 } 1352 } 1353 } 1354 } 1355 PetscCheck(total == Nbpt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrectly counted P trimmed polynomials"); 1356 PetscCall(PetscFree(deriv)); 1357 PetscCall(PetscFree(pattern)); 1358 PetscCall(PetscFree(form_atoms)); 1359 PetscCall(PetscFree(p_scalar)); 1360 PetscFunctionReturn(PETSC_SUCCESS); 1361 } 1362 1363 /*@ 1364 PetscDTPTrimmedEvalJet - Evaluate the jet (function and derivatives) of a basis of the trimmed polynomial k-forms up to 1365 a given degree. 1366 1367 Input Parameters: 1368 + dim - the number of variables in the multivariate polynomials 1369 . npoints - the number of points to evaluate the polynomials at 1370 . points - [npoints x dim] array of point coordinates 1371 . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space to evaluate. 1372 There are ((dim + degree) choose (dim + formDegree)) x ((degree + formDegree - 1) choose (formDegree)) polynomials in this space. 1373 (You can use `PetscDTPTrimmedSize()` to compute this size.) 1374 . formDegree - the degree of the form 1375 - jetDegree - the maximum order partial derivative to evaluate in the jet. There are ((dim + jetDegree) choose dim) partial derivatives 1376 in the jet. Choosing jetDegree = 0 means to evaluate just the function and no derivatives 1377 1378 Output Parameter: 1379 . p - an array containing the evaluations of the PKD polynomials' jets on the points. 1380 1381 Level: advanced 1382 1383 Notes: 1384 The size of `p` is `PetscDTPTrimmedSize()` x ((dim + formDegree) choose dim) x ((dim + k) 1385 choose dim) x npoints,which also describes the order of the dimensions of this 1386 four-dimensional array\: 1387 1388 the first (slowest varying) dimension is basis function index; 1389 the second dimension is component of the form; 1390 the third dimension is jet index; 1391 the fourth (fastest varying) dimension is the index of the evaluation point. 1392 1393 The ordering of the basis functions is not graded, so the basis functions are not nested by degree like `PetscDTPKDEvalJet()`. 1394 The basis functions are not an L2-orthonormal basis on any particular domain. 1395 1396 The implementation is based on the description of the trimmed polynomials up to degree r as 1397 the direct sum of polynomials up to degree (r-1) and the Koszul differential applied to 1398 homogeneous polynomials of degree (r-1). 1399 1400 .seealso: `PetscDTPKDEvalJet()`, `PetscDTPTrimmedSize()` 1401 @*/ 1402 PetscErrorCode PetscDTPTrimmedEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[]) 1403 { 1404 PetscFunctionBegin; 1405 PetscCall(PetscDTPTrimmedEvalJet_Internal(dim, npoints, points, degree, formDegree, jetDegree, p)); 1406 PetscFunctionReturn(PETSC_SUCCESS); 1407 } 1408 1409 /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V 1410 * with lds n; diag and subdiag are overwritten */ 1411 static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[], PetscReal eigs[], PetscScalar V[]) 1412 { 1413 char jobz = 'V'; /* eigenvalues and eigenvectors */ 1414 char range = 'A'; /* all eigenvalues will be found */ 1415 PetscReal VL = 0.; /* ignored because range is 'A' */ 1416 PetscReal VU = 0.; /* ignored because range is 'A' */ 1417 PetscBLASInt IL = 0; /* ignored because range is 'A' */ 1418 PetscBLASInt IU = 0; /* ignored because range is 'A' */ 1419 PetscReal abstol = 0.; /* unused */ 1420 PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */ 1421 PetscBLASInt *isuppz; 1422 PetscBLASInt lwork, liwork; 1423 PetscReal workquery; 1424 PetscBLASInt iworkquery; 1425 PetscBLASInt *iwork; 1426 PetscBLASInt info; 1427 PetscReal *work = NULL; 1428 1429 PetscFunctionBegin; 1430 #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG) 1431 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found"); 1432 #endif 1433 PetscCall(PetscBLASIntCast(n, &bn)); 1434 PetscCall(PetscBLASIntCast(n, &ldz)); 1435 #if !defined(PETSC_MISSING_LAPACK_STEGR) 1436 PetscCall(PetscMalloc1(2 * n, &isuppz)); 1437 lwork = -1; 1438 liwork = -1; 1439 PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, &workquery, &lwork, &iworkquery, &liwork, &info)); 1440 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEGR error"); 1441 lwork = (PetscBLASInt)workquery; 1442 liwork = (PetscBLASInt)iworkquery; 1443 PetscCall(PetscMalloc2(lwork, &work, liwork, &iwork)); 1444 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 1445 PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, work, &lwork, iwork, &liwork, &info)); 1446 PetscCall(PetscFPTrapPop()); 1447 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEGR error"); 1448 PetscCall(PetscFree2(work, iwork)); 1449 PetscCall(PetscFree(isuppz)); 1450 #elif !defined(PETSC_MISSING_LAPACK_STEQR) 1451 jobz = 'I'; /* Compute eigenvalues and eigenvectors of the 1452 tridiagonal matrix. Z is initialized to the identity 1453 matrix. */ 1454 PetscCall(PetscMalloc1(PetscMax(1, 2 * n - 2), &work)); 1455 PetscCallBLAS("LAPACKsteqr", LAPACKsteqr_("I", &bn, diag, subdiag, V, &ldz, work, &info)); 1456 PetscCall(PetscFPTrapPop()); 1457 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEQR error"); 1458 PetscCall(PetscFree(work)); 1459 PetscCall(PetscArraycpy(eigs, diag, n)); 1460 #endif 1461 PetscFunctionReturn(PETSC_SUCCESS); 1462 } 1463 1464 /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi 1465 * quadrature rules on the interval [-1, 1] */ 1466 static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw) 1467 { 1468 PetscReal twoab1; 1469 PetscInt m = n - 2; 1470 PetscReal a = alpha + 1.; 1471 PetscReal b = beta + 1.; 1472 PetscReal gra, grb; 1473 1474 PetscFunctionBegin; 1475 twoab1 = PetscPowReal(2., a + b - 1.); 1476 #if defined(PETSC_HAVE_LGAMMA) 1477 grb = PetscExpReal(2. * PetscLGamma(b + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + a + 1.) - (PetscLGamma(m + b + 1) + PetscLGamma(m + a + b + 1.))); 1478 gra = PetscExpReal(2. * PetscLGamma(a + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + b + 1.) - (PetscLGamma(m + a + 1) + PetscLGamma(m + a + b + 1.))); 1479 #else 1480 { 1481 PetscInt alphai = (PetscInt)alpha; 1482 PetscInt betai = (PetscInt)beta; 1483 1484 if ((PetscReal)alphai == alpha && (PetscReal)betai == beta) { 1485 PetscReal binom1, binom2; 1486 1487 PetscCall(PetscDTBinomial(m + b, b, &binom1)); 1488 PetscCall(PetscDTBinomial(m + a + b, b, &binom2)); 1489 grb = 1. / (binom1 * binom2); 1490 PetscCall(PetscDTBinomial(m + a, a, &binom1)); 1491 PetscCall(PetscDTBinomial(m + a + b, a, &binom2)); 1492 gra = 1. / (binom1 * binom2); 1493 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable."); 1494 } 1495 #endif 1496 *leftw = twoab1 * grb / b; 1497 *rightw = twoab1 * gra / a; 1498 PetscFunctionReturn(PETSC_SUCCESS); 1499 } 1500 1501 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x. 1502 Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */ 1503 static inline PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 1504 { 1505 PetscReal pn1, pn2; 1506 PetscReal cnm1, cnm1x, cnm2; 1507 PetscInt k; 1508 1509 PetscFunctionBegin; 1510 if (!n) { 1511 *P = 1.0; 1512 PetscFunctionReturn(PETSC_SUCCESS); 1513 } 1514 PetscDTJacobiRecurrence_Internal(1, a, b, cnm1, cnm1x, cnm2); 1515 pn2 = 1.; 1516 pn1 = cnm1 + cnm1x * x; 1517 if (n == 1) { 1518 *P = pn1; 1519 PetscFunctionReturn(PETSC_SUCCESS); 1520 } 1521 *P = 0.0; 1522 for (k = 2; k < n + 1; ++k) { 1523 PetscDTJacobiRecurrence_Internal(k, a, b, cnm1, cnm1x, cnm2); 1524 1525 *P = (cnm1 + cnm1x * x) * pn1 - cnm2 * pn2; 1526 pn2 = pn1; 1527 pn1 = *P; 1528 } 1529 PetscFunctionReturn(PETSC_SUCCESS); 1530 } 1531 1532 /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */ 1533 static inline PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P) 1534 { 1535 PetscReal nP; 1536 PetscInt i; 1537 1538 PetscFunctionBegin; 1539 *P = 0.0; 1540 if (k > n) PetscFunctionReturn(PETSC_SUCCESS); 1541 PetscCall(PetscDTComputeJacobi(a + k, b + k, n - k, x, &nP)); 1542 for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5; 1543 *P = nP; 1544 PetscFunctionReturn(PETSC_SUCCESS); 1545 } 1546 1547 static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[]) 1548 { 1549 PetscInt maxIter = 100; 1550 PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON)); 1551 PetscReal a1, a6, gf; 1552 PetscInt k; 1553 1554 PetscFunctionBegin; 1555 1556 a1 = PetscPowReal(2.0, a + b + 1); 1557 #if defined(PETSC_HAVE_LGAMMA) 1558 { 1559 PetscReal a2, a3, a4, a5; 1560 a2 = PetscLGamma(a + npoints + 1); 1561 a3 = PetscLGamma(b + npoints + 1); 1562 a4 = PetscLGamma(a + b + npoints + 1); 1563 a5 = PetscLGamma(npoints + 1); 1564 gf = PetscExpReal(a2 + a3 - (a4 + a5)); 1565 } 1566 #else 1567 { 1568 PetscInt ia, ib; 1569 1570 ia = (PetscInt)a; 1571 ib = (PetscInt)b; 1572 gf = 1.; 1573 if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */ 1574 for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k); 1575 } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */ 1576 for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k); 1577 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable."); 1578 } 1579 #endif 1580 1581 a6 = a1 * gf; 1582 /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses. 1583 Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */ 1584 for (k = 0; k < npoints; ++k) { 1585 PetscReal r = PetscCosReal(PETSC_PI * (1. - (4. * k + 3. + 2. * b) / (4. * npoints + 2. * (a + b + 1.)))), dP; 1586 PetscInt j; 1587 1588 if (k > 0) r = 0.5 * (r + x[k - 1]); 1589 for (j = 0; j < maxIter; ++j) { 1590 PetscReal s = 0.0, delta, f, fp; 1591 PetscInt i; 1592 1593 for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]); 1594 PetscCall(PetscDTComputeJacobi(a, b, npoints, r, &f)); 1595 PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp)); 1596 delta = f / (fp - f * s); 1597 r = r - delta; 1598 if (PetscAbsReal(delta) < eps) break; 1599 } 1600 x[k] = r; 1601 PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP)); 1602 w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP); 1603 } 1604 PetscFunctionReturn(PETSC_SUCCESS); 1605 } 1606 1607 /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi 1608 * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */ 1609 static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s) 1610 { 1611 PetscInt i; 1612 1613 PetscFunctionBegin; 1614 for (i = 0; i < nPoints; i++) { 1615 PetscReal A, B, C; 1616 1617 PetscDTJacobiRecurrence_Internal(i + 1, a, b, A, B, C); 1618 d[i] = -A / B; 1619 if (i) s[i - 1] *= C / B; 1620 if (i < nPoints - 1) s[i] = 1. / B; 1621 } 1622 PetscFunctionReturn(PETSC_SUCCESS); 1623 } 1624 1625 static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w) 1626 { 1627 PetscReal mu0; 1628 PetscReal ga, gb, gab; 1629 PetscInt i; 1630 1631 PetscFunctionBegin; 1632 PetscCall(PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite)); 1633 1634 #if defined(PETSC_HAVE_TGAMMA) 1635 ga = PetscTGamma(a + 1); 1636 gb = PetscTGamma(b + 1); 1637 gab = PetscTGamma(a + b + 2); 1638 #else 1639 { 1640 PetscInt ia, ib; 1641 1642 ia = (PetscInt)a; 1643 ib = (PetscInt)b; 1644 if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */ 1645 PetscCall(PetscDTFactorial(ia, &ga)); 1646 PetscCall(PetscDTFactorial(ib, &gb)); 1647 PetscCall(PetscDTFactorial(ia + ib + 1, &gab)); 1648 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "tgamma() - math routine is unavailable."); 1649 } 1650 #endif 1651 mu0 = PetscPowReal(2., a + b + 1.) * ga * gb / gab; 1652 1653 #if defined(PETSCDTGAUSSIANQUADRATURE_EIG) 1654 { 1655 PetscReal *diag, *subdiag; 1656 PetscScalar *V; 1657 1658 PetscCall(PetscMalloc2(npoints, &diag, npoints, &subdiag)); 1659 PetscCall(PetscMalloc1(npoints * npoints, &V)); 1660 PetscCall(PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag)); 1661 for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]); 1662 PetscCall(PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V)); 1663 for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0; 1664 PetscCall(PetscFree(V)); 1665 PetscCall(PetscFree2(diag, subdiag)); 1666 } 1667 #else 1668 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found"); 1669 #endif 1670 { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the 1671 eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that 1672 the eigenvalues are sorted */ 1673 PetscBool sorted; 1674 1675 PetscCall(PetscSortedReal(npoints, x, &sorted)); 1676 if (!sorted) { 1677 PetscInt *order, i; 1678 PetscReal *tmp; 1679 1680 PetscCall(PetscMalloc2(npoints, &order, npoints, &tmp)); 1681 for (i = 0; i < npoints; i++) order[i] = i; 1682 PetscCall(PetscSortRealWithPermutation(npoints, x, order)); 1683 PetscCall(PetscArraycpy(tmp, x, npoints)); 1684 for (i = 0; i < npoints; i++) x[i] = tmp[order[i]]; 1685 PetscCall(PetscArraycpy(tmp, w, npoints)); 1686 for (i = 0; i < npoints; i++) w[i] = tmp[order[i]]; 1687 PetscCall(PetscFree2(order, tmp)); 1688 } 1689 } 1690 PetscFunctionReturn(PETSC_SUCCESS); 1691 } 1692 1693 static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1694 { 1695 PetscFunctionBegin; 1696 PetscCheck(npoints >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive"); 1697 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1698 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1."); 1699 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1."); 1700 1701 if (newton) PetscCall(PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w)); 1702 else PetscCall(PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w)); 1703 if (alpha == beta) { /* symmetrize */ 1704 PetscInt i; 1705 for (i = 0; i < (npoints + 1) / 2; i++) { 1706 PetscInt j = npoints - 1 - i; 1707 PetscReal xi = x[i]; 1708 PetscReal xj = x[j]; 1709 PetscReal wi = w[i]; 1710 PetscReal wj = w[j]; 1711 1712 x[i] = (xi - xj) / 2.; 1713 x[j] = (xj - xi) / 2.; 1714 w[i] = w[j] = (wi + wj) / 2.; 1715 } 1716 } 1717 PetscFunctionReturn(PETSC_SUCCESS); 1718 } 1719 1720 /*@ 1721 PetscDTGaussJacobiQuadrature - quadrature for the interval [a, b] with the weight function 1722 $(x-a)^\alpha (x-b)^\beta$. 1723 1724 Not Collective 1725 1726 Input Parameters: 1727 + npoints - the number of points in the quadrature rule 1728 . a - the left endpoint of the interval 1729 . b - the right endpoint of the interval 1730 . alpha - the left exponent 1731 - beta - the right exponent 1732 1733 Output Parameters: 1734 + x - array of length `npoints`, the locations of the quadrature points 1735 - w - array of length `npoints`, the weights of the quadrature points 1736 1737 Level: intermediate 1738 1739 Note: 1740 This quadrature rule is exact for polynomials up to degree 2*npoints - 1. 1741 1742 .seealso: `PetscDTGaussQuadrature()` 1743 @*/ 1744 PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1745 { 1746 PetscInt i; 1747 1748 PetscFunctionBegin; 1749 PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal)); 1750 if (a != -1. || b != 1.) { /* shift */ 1751 for (i = 0; i < npoints; i++) { 1752 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1753 w[i] *= (b - a) / 2.; 1754 } 1755 } 1756 PetscFunctionReturn(PETSC_SUCCESS); 1757 } 1758 1759 static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1760 { 1761 PetscInt i; 1762 1763 PetscFunctionBegin; 1764 PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive"); 1765 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1766 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1."); 1767 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1."); 1768 1769 x[0] = -1.; 1770 x[npoints - 1] = 1.; 1771 if (npoints > 2) PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints - 2, alpha + 1., beta + 1., &x[1], &w[1], newton)); 1772 for (i = 1; i < npoints - 1; i++) w[i] /= (1. - x[i] * x[i]); 1773 PetscCall(PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints - 1])); 1774 PetscFunctionReturn(PETSC_SUCCESS); 1775 } 1776 1777 /*@ 1778 PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval [a, b] with the weight function 1779 $(x-a)^\alpha (x-b)^\beta$, with endpoints a and b included as quadrature points. 1780 1781 Not Collective 1782 1783 Input Parameters: 1784 + npoints - the number of points in the quadrature rule 1785 . a - the left endpoint of the interval 1786 . b - the right endpoint of the interval 1787 . alpha - the left exponent 1788 - beta - the right exponent 1789 1790 Output Parameters: 1791 + x - array of length `npoints`, the locations of the quadrature points 1792 - w - array of length `npoints`, the weights of the quadrature points 1793 1794 Level: intermediate 1795 1796 Note: 1797 This quadrature rule is exact for polynomials up to degree 2*npoints - 3. 1798 1799 .seealso: `PetscDTGaussJacobiQuadrature()` 1800 @*/ 1801 PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1802 { 1803 PetscInt i; 1804 1805 PetscFunctionBegin; 1806 PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal)); 1807 if (a != -1. || b != 1.) { /* shift */ 1808 for (i = 0; i < npoints; i++) { 1809 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1810 w[i] *= (b - a) / 2.; 1811 } 1812 } 1813 PetscFunctionReturn(PETSC_SUCCESS); 1814 } 1815 1816 /*@ 1817 PetscDTGaussQuadrature - create Gauss-Legendre quadrature 1818 1819 Not Collective 1820 1821 Input Parameters: 1822 + npoints - number of points 1823 . a - left end of interval (often-1) 1824 - b - right end of interval (often +1) 1825 1826 Output Parameters: 1827 + x - quadrature points 1828 - w - quadrature weights 1829 1830 Level: intermediate 1831 1832 References: 1833 . * - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969. 1834 1835 .seealso: `PetscDTLegendreEval()`, `PetscDTGaussJacobiQuadrature()` 1836 @*/ 1837 PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w) 1838 { 1839 PetscInt i; 1840 1841 PetscFunctionBegin; 1842 PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal)); 1843 if (a != -1. || b != 1.) { /* shift */ 1844 for (i = 0; i < npoints; i++) { 1845 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1846 w[i] *= (b - a) / 2.; 1847 } 1848 } 1849 PetscFunctionReturn(PETSC_SUCCESS); 1850 } 1851 1852 /*@C 1853 PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre 1854 nodes of a given size on the domain [-1,1] 1855 1856 Not Collective 1857 1858 Input Parameters: 1859 + npoints - number of grid nodes 1860 - type - `PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA` or `PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON` 1861 1862 Output Parameters: 1863 + x - quadrature points 1864 - w - quadrature weights 1865 1866 Level: intermediate 1867 1868 Notes: 1869 For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not 1870 close enough to the desired solution 1871 1872 These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes 1873 1874 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 1875 1876 .seealso: `PetscDTGaussQuadrature()`, `PetscGaussLobattoLegendreCreateType` 1877 1878 @*/ 1879 PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints, PetscGaussLobattoLegendreCreateType type, PetscReal *x, PetscReal *w) 1880 { 1881 PetscBool newton; 1882 1883 PetscFunctionBegin; 1884 PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide at least 2 grid points per element"); 1885 newton = (PetscBool)(type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON); 1886 PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton)); 1887 PetscFunctionReturn(PETSC_SUCCESS); 1888 } 1889 1890 /*@ 1891 PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature 1892 1893 Not Collective 1894 1895 Input Parameters: 1896 + dim - The spatial dimension 1897 . Nc - The number of components 1898 . npoints - number of points in one dimension 1899 . a - left end of interval (often-1) 1900 - b - right end of interval (often +1) 1901 1902 Output Parameter: 1903 . q - A `PetscQuadrature` object 1904 1905 Level: intermediate 1906 1907 .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()` 1908 @*/ 1909 PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 1910 { 1911 DMPolytopeType ct; 1912 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints * PetscSqr(npoints) : PetscSqr(npoints) : npoints; 1913 PetscReal *x, *w, *xw, *ww; 1914 1915 PetscFunctionBegin; 1916 PetscCall(PetscMalloc1(totpoints * dim, &x)); 1917 PetscCall(PetscMalloc1(totpoints * Nc, &w)); 1918 /* Set up the Golub-Welsch system */ 1919 switch (dim) { 1920 case 0: 1921 ct = DM_POLYTOPE_POINT; 1922 PetscCall(PetscFree(x)); 1923 PetscCall(PetscFree(w)); 1924 PetscCall(PetscMalloc1(1, &x)); 1925 PetscCall(PetscMalloc1(Nc, &w)); 1926 totpoints = 1; 1927 x[0] = 0.0; 1928 for (PetscInt c = 0; c < Nc; ++c) w[c] = 1.0; 1929 break; 1930 case 1: 1931 ct = DM_POLYTOPE_SEGMENT; 1932 PetscCall(PetscMalloc1(npoints, &ww)); 1933 PetscCall(PetscDTGaussQuadrature(npoints, a, b, x, ww)); 1934 for (PetscInt i = 0; i < npoints; ++i) 1935 for (PetscInt c = 0; c < Nc; ++c) w[i * Nc + c] = ww[i]; 1936 PetscCall(PetscFree(ww)); 1937 break; 1938 case 2: 1939 ct = DM_POLYTOPE_QUADRILATERAL; 1940 PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww)); 1941 PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww)); 1942 for (PetscInt i = 0; i < npoints; ++i) { 1943 for (PetscInt j = 0; j < npoints; ++j) { 1944 x[(i * npoints + j) * dim + 0] = xw[i]; 1945 x[(i * npoints + j) * dim + 1] = xw[j]; 1946 for (PetscInt c = 0; c < Nc; ++c) w[(i * npoints + j) * Nc + c] = ww[i] * ww[j]; 1947 } 1948 } 1949 PetscCall(PetscFree2(xw, ww)); 1950 break; 1951 case 3: 1952 ct = DM_POLYTOPE_HEXAHEDRON; 1953 PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww)); 1954 PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww)); 1955 for (PetscInt i = 0; i < npoints; ++i) { 1956 for (PetscInt j = 0; j < npoints; ++j) { 1957 for (PetscInt k = 0; k < npoints; ++k) { 1958 x[((i * npoints + j) * npoints + k) * dim + 0] = xw[i]; 1959 x[((i * npoints + j) * npoints + k) * dim + 1] = xw[j]; 1960 x[((i * npoints + j) * npoints + k) * dim + 2] = xw[k]; 1961 for (PetscInt c = 0; c < Nc; ++c) w[((i * npoints + j) * npoints + k) * Nc + c] = ww[i] * ww[j] * ww[k]; 1962 } 1963 } 1964 } 1965 PetscCall(PetscFree2(xw, ww)); 1966 break; 1967 default: 1968 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %" PetscInt_FMT, dim); 1969 } 1970 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 1971 PetscCall(PetscQuadratureSetCellType(*q, ct)); 1972 PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1)); 1973 PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w)); 1974 PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "GaussTensor")); 1975 PetscFunctionReturn(PETSC_SUCCESS); 1976 } 1977 1978 /*@ 1979 PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex 1980 1981 Not Collective 1982 1983 Input Parameters: 1984 + dim - The simplex dimension 1985 . Nc - The number of components 1986 . npoints - The number of points in one dimension 1987 . a - left end of interval (often-1) 1988 - b - right end of interval (often +1) 1989 1990 Output Parameter: 1991 . q - A `PetscQuadrature` object 1992 1993 Level: intermediate 1994 1995 Note: 1996 For `dim` == 1, this is Gauss-Legendre quadrature 1997 1998 References: 1999 . * - Karniadakis and Sherwin. FIAT 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 /*@C 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 = DMPolytopeTypeGetNumArrangments(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