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 PetscValidPointer(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 PetscValidPointer(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 PetscValidPointer(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 PetscValidIntPointer(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 PetscValidIntPointer(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 Note: 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 PetscValidIntPointer(dim, 2); 310 *dim = q->dim; 311 } 312 if (Nc) { 313 PetscValidIntPointer(Nc, 3); 314 *Nc = q->Nc; 315 } 316 if (npoints) { 317 PetscValidIntPointer(npoints, 4); 318 *npoints = q->numPoints; 319 } 320 if (points) { 321 PetscValidPointer(points, 5); 322 *points = q->points; 323 } 324 if (weights) { 325 PetscValidPointer(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 PetscValidBoolPointer(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 PetscValidRealPointer(points, 5); 569 q->points = points; 570 } 571 if (weights) { 572 PetscValidRealPointer(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 . dim - 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 PetscValidRealPointer(v0, 3); 668 PetscValidRealPointer(jac, 4); 669 PetscValidPointer(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 The Jacobi polynomials with indices $\alpha$ and $\beta$ are orthogonal with respect to the weighted inner product 837 $\langle f, g \rangle = \int_{-1}^1 (1+x)^{\alpha} (1-x)^{\beta} f(x) g(x) dx$. 838 839 Input Parameters: 840 + alpha - the left exponent of the weight 841 . beta - the right exponetn of the weight 842 . npoints - the number of points to evaluate the polynomials at 843 . points - [npoints] array of point coordinates 844 . degree - the maximm degree polynomial space to evaluate, (degree + 1) will be evaluated total. 845 - k - the maximum derivative to evaluate in the jet, (k + 1) will be evaluated total. 846 847 Output Parameter: 848 . p - an array containing the evaluations of the Jacobi polynomials's jets on the points. the size is (degree + 1) x 849 (k + 1) x npoints, which also describes the order of the dimensions of this three-dimensional array: the first 850 (slowest varying) dimension is polynomial degree; the second dimension is derivative order; the third (fastest 851 varying) dimension is the index of the evaluation point. 852 853 Level: advanced 854 855 .seealso: `PetscDTJacobiEval()`, `PetscDTPKDEvalJet()` 856 @*/ 857 PetscErrorCode PetscDTJacobiEvalJet(PetscReal alpha, PetscReal beta, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[]) 858 { 859 PetscInt i, j, l; 860 PetscInt *degrees; 861 PetscReal *psingle; 862 863 PetscFunctionBegin; 864 if (degree == 0) { 865 PetscInt zero = 0; 866 867 for (i = 0; i <= k; i++) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, 1, &zero, &p[i * npoints])); 868 PetscFunctionReturn(PETSC_SUCCESS); 869 } 870 PetscCall(PetscMalloc1(degree + 1, °rees)); 871 PetscCall(PetscMalloc1((degree + 1) * npoints, &psingle)); 872 for (i = 0; i <= degree; i++) degrees[i] = i; 873 for (i = 0; i <= k; i++) { 874 PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, degree + 1, degrees, psingle)); 875 for (j = 0; j <= degree; j++) { 876 for (l = 0; l < npoints; l++) p[(j * (k + 1) + i) * npoints + l] = psingle[l * (degree + 1) + j]; 877 } 878 } 879 PetscCall(PetscFree(psingle)); 880 PetscCall(PetscFree(degrees)); 881 PetscFunctionReturn(PETSC_SUCCESS); 882 } 883 884 /*@ 885 PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$ at a set of points 886 at points 887 888 Not Collective 889 890 Input Parameters: 891 + npoints - number of spatial points to evaluate at 892 . alpha - the left exponent > -1 893 . beta - the right exponent > -1 894 . points - array of locations to evaluate at 895 . ndegree - number of basis degrees to evaluate 896 - degrees - sorted array of degrees to evaluate 897 898 Output Parameters: 899 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL) 900 . D - row-oriented derivative evaluation matrix (or NULL) 901 - D2 - row-oriented second derivative evaluation matrix (or NULL) 902 903 Level: intermediate 904 905 .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()` 906 @*/ 907 PetscErrorCode PetscDTJacobiEval(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *B, PetscReal *D, PetscReal *D2) 908 { 909 PetscFunctionBegin; 910 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1."); 911 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1."); 912 if (!npoints || !ndegree) PetscFunctionReturn(PETSC_SUCCESS); 913 if (B) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B)); 914 if (D) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D)); 915 if (D2) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2)); 916 PetscFunctionReturn(PETSC_SUCCESS); 917 } 918 919 /*@ 920 PetscDTLegendreEval - evaluate Legendre polynomials at points 921 922 Not Collective 923 924 Input Parameters: 925 + npoints - number of spatial points to evaluate at 926 . points - array of locations to evaluate at 927 . ndegree - number of basis degrees to evaluate 928 - degrees - sorted array of degrees to evaluate 929 930 Output Parameters: 931 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL) 932 . D - row-oriented derivative evaluation matrix (or NULL) 933 - D2 - row-oriented second derivative evaluation matrix (or NULL) 934 935 Level: intermediate 936 937 .seealso: `PetscDTGaussQuadrature()` 938 @*/ 939 PetscErrorCode PetscDTLegendreEval(PetscInt npoints, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *B, PetscReal *D, PetscReal *D2) 940 { 941 PetscFunctionBegin; 942 PetscCall(PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2)); 943 PetscFunctionReturn(PETSC_SUCCESS); 944 } 945 946 /*@ 947 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) 948 949 Input Parameters: 950 + len - the desired length of the degree tuple 951 - index - the index to convert: should be >= 0 952 953 Output Parameter: 954 . degtup - will be filled with a tuple of degrees 955 956 Level: beginner 957 958 Note: 959 For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples 960 acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the 961 last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1). 962 963 .seealso: `PetscDTGradedOrderToIndex()` 964 @*/ 965 PetscErrorCode PetscDTIndexToGradedOrder(PetscInt len, PetscInt index, PetscInt degtup[]) 966 { 967 PetscInt i, total; 968 PetscInt sum; 969 970 PetscFunctionBeginHot; 971 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 972 PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative"); 973 total = 1; 974 sum = 0; 975 while (index >= total) { 976 index -= total; 977 total = (total * (len + sum)) / (sum + 1); 978 sum++; 979 } 980 for (i = 0; i < len; i++) { 981 PetscInt c; 982 983 degtup[i] = sum; 984 for (c = 0, total = 1; c < sum; c++) { 985 /* going into the loop, total is the number of way to have a tuple of sum exactly c with length len - 1 - i */ 986 if (index < total) break; 987 index -= total; 988 total = (total * (len - 1 - i + c)) / (c + 1); 989 degtup[i]--; 990 } 991 sum -= degtup[i]; 992 } 993 PetscFunctionReturn(PETSC_SUCCESS); 994 } 995 996 /*@ 997 PetscDTGradedOrderToIndex - convert a tuple into an index in a graded order, the inverse of `PetscDTIndexToGradedOrder()`. 998 999 Input Parameters: 1000 + len - the length of the degree tuple 1001 - degtup - tuple with this length 1002 1003 Output Parameter: 1004 . index - index in graded order: >= 0 1005 1006 Level: Beginner 1007 1008 Note: 1009 For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples 1010 acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the 1011 last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1). 1012 1013 .seealso: `PetscDTIndexToGradedOrder()` 1014 @*/ 1015 PetscErrorCode PetscDTGradedOrderToIndex(PetscInt len, const PetscInt degtup[], PetscInt *index) 1016 { 1017 PetscInt i, idx, sum, total; 1018 1019 PetscFunctionBeginHot; 1020 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 1021 for (i = 0, sum = 0; i < len; i++) sum += degtup[i]; 1022 idx = 0; 1023 total = 1; 1024 for (i = 0; i < sum; i++) { 1025 idx += total; 1026 total = (total * (len + i)) / (i + 1); 1027 } 1028 for (i = 0; i < len - 1; i++) { 1029 PetscInt c; 1030 1031 total = 1; 1032 sum -= degtup[i]; 1033 for (c = 0; c < sum; c++) { 1034 idx += total; 1035 total = (total * (len - 1 - i + c)) / (c + 1); 1036 } 1037 } 1038 *index = idx; 1039 PetscFunctionReturn(PETSC_SUCCESS); 1040 } 1041 1042 static PetscBool PKDCite = PETSC_FALSE; 1043 const char PKDCitation[] = "@article{Kirby2010,\n" 1044 " title={Singularity-free evaluation of collapsed-coordinate orthogonal polynomials},\n" 1045 " author={Kirby, Robert C},\n" 1046 " journal={ACM Transactions on Mathematical Software (TOMS)},\n" 1047 " volume={37},\n" 1048 " number={1},\n" 1049 " pages={1--16},\n" 1050 " year={2010},\n" 1051 " publisher={ACM New York, NY, USA}\n}\n"; 1052 1053 /*@ 1054 PetscDTPKDEvalJet - Evaluate the jet (function and derivatives) of the Proriol-Koornwinder-Dubiner (PKD) basis for 1055 the space of polynomials up to a given degree. The PKD basis is L2-orthonormal on the biunit simplex (which is used 1056 as the reference element for finite elements in PETSc), which makes it a stable basis to use for evaluating 1057 polynomials in that domain. 1058 1059 Input Parameters: 1060 + dim - the number of variables in the multivariate polynomials 1061 . npoints - the number of points to evaluate the polynomials at 1062 . points - [npoints x dim] array of point coordinates 1063 . 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. 1064 - k - the maximum order partial derivative to evaluate in the jet. There are (dim + k choose dim) partial derivatives 1065 in the jet. Choosing k = 0 means to evaluate just the function and no derivatives 1066 1067 Output Parameter: 1068 . p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is ((dim + degree) 1069 choose dim) x ((dim + k) choose dim) x npoints, which also describes the order of the dimensions of this 1070 three-dimensional array: the first (slowest varying) dimension is basis function index; the second dimension is jet 1071 index; the third (fastest varying) dimension is the index of the evaluation point. 1072 1073 Level: advanced 1074 1075 Notes: 1076 The ordering of the basis functions, and the ordering of the derivatives in the jet, both follow the graded 1077 ordering of `PetscDTIndexToGradedOrder()` and `PetscDTGradedOrderToIndex()`. For example, in 3D, the polynomial with 1078 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); 1079 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). 1080 1081 The implementation uses Kirby's singularity-free evaluation algorithm, https://doi.org/10.1145/1644001.1644006. 1082 1083 .seealso: `PetscDTGradedOrderToIndex()`, `PetscDTIndexToGradedOrder()`, `PetscDTJacobiEvalJet()` 1084 @*/ 1085 PetscErrorCode PetscDTPKDEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[]) 1086 { 1087 PetscInt degidx, kidx, d, pt; 1088 PetscInt Nk, Ndeg; 1089 PetscInt *ktup, *degtup; 1090 PetscReal *scales, initscale, scaleexp; 1091 1092 PetscFunctionBegin; 1093 PetscCall(PetscCitationsRegister(PKDCitation, &PKDCite)); 1094 PetscCall(PetscDTBinomialInt(dim + k, k, &Nk)); 1095 PetscCall(PetscDTBinomialInt(degree + dim, degree, &Ndeg)); 1096 PetscCall(PetscMalloc2(dim, °tup, dim, &ktup)); 1097 PetscCall(PetscMalloc1(Ndeg, &scales)); 1098 initscale = 1.; 1099 if (dim > 1) { 1100 PetscCall(PetscDTBinomial(dim, 2, &scaleexp)); 1101 initscale = PetscPowReal(2., scaleexp * 0.5); 1102 } 1103 for (degidx = 0; degidx < Ndeg; degidx++) { 1104 PetscInt e, i; 1105 PetscInt m1idx = -1, m2idx = -1; 1106 PetscInt n; 1107 PetscInt degsum; 1108 PetscReal alpha; 1109 PetscReal cnm1, cnm1x, cnm2; 1110 PetscReal norm; 1111 1112 PetscCall(PetscDTIndexToGradedOrder(dim, degidx, degtup)); 1113 for (d = dim - 1; d >= 0; d--) 1114 if (degtup[d]) break; 1115 if (d < 0) { /* constant is 1 everywhere, all derivatives are zero */ 1116 scales[degidx] = initscale; 1117 for (e = 0; e < dim; e++) { 1118 PetscCall(PetscDTJacobiNorm(e, 0., 0, &norm)); 1119 scales[degidx] /= norm; 1120 } 1121 for (i = 0; i < npoints; i++) p[degidx * Nk * npoints + i] = 1.; 1122 for (i = 0; i < (Nk - 1) * npoints; i++) p[(degidx * Nk + 1) * npoints + i] = 0.; 1123 continue; 1124 } 1125 n = degtup[d]; 1126 degtup[d]--; 1127 PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m1idx)); 1128 if (degtup[d] > 0) { 1129 degtup[d]--; 1130 PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m2idx)); 1131 degtup[d]++; 1132 } 1133 degtup[d]++; 1134 for (e = 0, degsum = 0; e < d; e++) degsum += degtup[e]; 1135 alpha = 2 * degsum + d; 1136 PetscDTJacobiRecurrence_Internal(n, alpha, 0., cnm1, cnm1x, cnm2); 1137 1138 scales[degidx] = initscale; 1139 for (e = 0, degsum = 0; e < dim; e++) { 1140 PetscInt f; 1141 PetscReal ealpha; 1142 PetscReal enorm; 1143 1144 ealpha = 2 * degsum + e; 1145 for (f = 0; f < degsum; f++) scales[degidx] *= 2.; 1146 PetscCall(PetscDTJacobiNorm(ealpha, 0., degtup[e], &enorm)); 1147 scales[degidx] /= enorm; 1148 degsum += degtup[e]; 1149 } 1150 1151 for (pt = 0; pt < npoints; pt++) { 1152 /* compute the multipliers */ 1153 PetscReal thetanm1, thetanm1x, thetanm2; 1154 1155 thetanm1x = dim - (d + 1) + 2. * points[pt * dim + d]; 1156 for (e = d + 1; e < dim; e++) thetanm1x += points[pt * dim + e]; 1157 thetanm1x *= 0.5; 1158 thetanm1 = (2. - (dim - (d + 1))); 1159 for (e = d + 1; e < dim; e++) thetanm1 -= points[pt * dim + e]; 1160 thetanm1 *= 0.5; 1161 thetanm2 = thetanm1 * thetanm1; 1162 1163 for (kidx = 0; kidx < Nk; kidx++) { 1164 PetscInt f; 1165 1166 PetscCall(PetscDTIndexToGradedOrder(dim, kidx, ktup)); 1167 /* first sum in the same derivative terms */ 1168 p[(degidx * Nk + kidx) * npoints + pt] = (cnm1 * thetanm1 + cnm1x * thetanm1x) * p[(m1idx * Nk + kidx) * npoints + pt]; 1169 if (m2idx >= 0) p[(degidx * Nk + kidx) * npoints + pt] -= cnm2 * thetanm2 * p[(m2idx * Nk + kidx) * npoints + pt]; 1170 1171 for (f = d; f < dim; f++) { 1172 PetscInt km1idx, mplty = ktup[f]; 1173 1174 if (!mplty) continue; 1175 ktup[f]--; 1176 PetscCall(PetscDTGradedOrderToIndex(dim, ktup, &km1idx)); 1177 1178 /* the derivative of cnm1x * thetanm1x wrt x variable f is 0.5 * cnm1x if f > d otherwise it is cnm1x */ 1179 /* the derivative of cnm1 * thetanm1 wrt x variable f is 0 if f == d, otherwise it is -0.5 * cnm1 */ 1180 /* the derivative of -cnm2 * thetanm2 wrt x variable f is 0 if f == d, otherwise it is cnm2 * thetanm1 */ 1181 if (f > d) { 1182 PetscInt f2; 1183 1184 p[(degidx * Nk + kidx) * npoints + pt] += mplty * 0.5 * (cnm1x - cnm1) * p[(m1idx * Nk + km1idx) * npoints + pt]; 1185 if (m2idx >= 0) { 1186 p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm2 * thetanm1 * p[(m2idx * Nk + km1idx) * npoints + pt]; 1187 /* second derivatives of -cnm2 * thetanm2 wrt x variable f,f2 is like - 0.5 * cnm2 */ 1188 for (f2 = f; f2 < dim; f2++) { 1189 PetscInt km2idx, mplty2 = ktup[f2]; 1190 PetscInt factor; 1191 1192 if (!mplty2) continue; 1193 ktup[f2]--; 1194 PetscCall(PetscDTGradedOrderToIndex(dim, ktup, &km2idx)); 1195 1196 factor = mplty * mplty2; 1197 if (f == f2) factor /= 2; 1198 p[(degidx * Nk + kidx) * npoints + pt] -= 0.5 * factor * cnm2 * p[(m2idx * Nk + km2idx) * npoints + pt]; 1199 ktup[f2]++; 1200 } 1201 } 1202 } else { 1203 p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm1x * p[(m1idx * Nk + km1idx) * npoints + pt]; 1204 } 1205 ktup[f]++; 1206 } 1207 } 1208 } 1209 } 1210 for (degidx = 0; degidx < Ndeg; degidx++) { 1211 PetscReal scale = scales[degidx]; 1212 PetscInt i; 1213 1214 for (i = 0; i < Nk * npoints; i++) p[degidx * Nk * npoints + i] *= scale; 1215 } 1216 PetscCall(PetscFree(scales)); 1217 PetscCall(PetscFree2(degtup, ktup)); 1218 PetscFunctionReturn(PETSC_SUCCESS); 1219 } 1220 1221 /*@ 1222 PetscDTPTrimmedSize - The size of the trimmed polynomial space of k-forms with a given degree and form degree, 1223 which can be evaluated in `PetscDTPTrimmedEvalJet()`. 1224 1225 Input Parameters: 1226 + dim - the number of variables in the multivariate polynomials 1227 . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space. 1228 - formDegree - the degree of the form 1229 1230 Output Parameter: 1231 - size - The number ((`dim` + `degree`) choose (`dim` + `formDegree`)) x ((`degree` + `formDegree` - 1) choose (`formDegree`)) 1232 1233 Level: advanced 1234 1235 .seealso: `PetscDTPTrimmedEvalJet()` 1236 @*/ 1237 PetscErrorCode PetscDTPTrimmedSize(PetscInt dim, PetscInt degree, PetscInt formDegree, PetscInt *size) 1238 { 1239 PetscInt Nrk, Nbpt; // number of trimmed polynomials 1240 1241 PetscFunctionBegin; 1242 formDegree = PetscAbsInt(formDegree); 1243 PetscCall(PetscDTBinomialInt(degree + dim, degree + formDegree, &Nbpt)); 1244 PetscCall(PetscDTBinomialInt(degree + formDegree - 1, formDegree, &Nrk)); 1245 Nbpt *= Nrk; 1246 *size = Nbpt; 1247 PetscFunctionReturn(PETSC_SUCCESS); 1248 } 1249 1250 /* there was a reference implementation based on section 4.4 of Arnold, Falk & Winther (acta numerica, 2006), but it 1251 * was inferior to this implementation */ 1252 static PetscErrorCode PetscDTPTrimmedEvalJet_Internal(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[]) 1253 { 1254 PetscInt formDegreeOrig = formDegree; 1255 PetscBool formNegative = (formDegreeOrig < 0) ? PETSC_TRUE : PETSC_FALSE; 1256 1257 PetscFunctionBegin; 1258 formDegree = PetscAbsInt(formDegreeOrig); 1259 if (formDegree == 0) { 1260 PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree, jetDegree, p)); 1261 PetscFunctionReturn(PETSC_SUCCESS); 1262 } 1263 if (formDegree == dim) { 1264 PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p)); 1265 PetscFunctionReturn(PETSC_SUCCESS); 1266 } 1267 PetscInt Nbpt; 1268 PetscCall(PetscDTPTrimmedSize(dim, degree, formDegree, &Nbpt)); 1269 PetscInt Nf; 1270 PetscCall(PetscDTBinomialInt(dim, formDegree, &Nf)); 1271 PetscInt Nk; 1272 PetscCall(PetscDTBinomialInt(dim + jetDegree, dim, &Nk)); 1273 PetscCall(PetscArrayzero(p, Nbpt * Nf * Nk * npoints)); 1274 1275 PetscInt Nbpm1; // number of scalar polynomials up to degree - 1; 1276 PetscCall(PetscDTBinomialInt(dim + degree - 1, dim, &Nbpm1)); 1277 PetscReal *p_scalar; 1278 PetscCall(PetscMalloc1(Nbpm1 * Nk * npoints, &p_scalar)); 1279 PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p_scalar)); 1280 PetscInt total = 0; 1281 // First add the full polynomials up to degree - 1 into the basis: take the scalar 1282 // and copy one for each form component 1283 for (PetscInt i = 0; i < Nbpm1; i++) { 1284 const PetscReal *src = &p_scalar[i * Nk * npoints]; 1285 for (PetscInt f = 0; f < Nf; f++) { 1286 PetscReal *dest = &p[(total++ * Nf + f) * Nk * npoints]; 1287 PetscCall(PetscArraycpy(dest, src, Nk * npoints)); 1288 } 1289 } 1290 PetscInt *form_atoms; 1291 PetscCall(PetscMalloc1(formDegree + 1, &form_atoms)); 1292 // construct the interior product pattern 1293 PetscInt(*pattern)[3]; 1294 PetscInt Nf1; // number of formDegree + 1 forms 1295 PetscCall(PetscDTBinomialInt(dim, formDegree + 1, &Nf1)); 1296 PetscInt nnz = Nf1 * (formDegree + 1); 1297 PetscCall(PetscMalloc1(Nf1 * (formDegree + 1), &pattern)); 1298 PetscCall(PetscDTAltVInteriorPattern(dim, formDegree + 1, pattern)); 1299 PetscReal centroid = (1. - dim) / (dim + 1.); 1300 PetscInt *deriv; 1301 PetscCall(PetscMalloc1(dim, &deriv)); 1302 for (PetscInt d = dim; d >= formDegree + 1; d--) { 1303 PetscInt Nfd1; // number of formDegree + 1 forms in dimension d that include dx_0 1304 // (equal to the number of formDegree forms in dimension d-1) 1305 PetscCall(PetscDTBinomialInt(d - 1, formDegree, &Nfd1)); 1306 // The number of homogeneous (degree-1) scalar polynomials in d variables 1307 PetscInt Nh; 1308 PetscCall(PetscDTBinomialInt(d - 1 + degree - 1, d - 1, &Nh)); 1309 const PetscReal *h_scalar = &p_scalar[(Nbpm1 - Nh) * Nk * npoints]; 1310 for (PetscInt b = 0; b < Nh; b++) { 1311 const PetscReal *h_s = &h_scalar[b * Nk * npoints]; 1312 for (PetscInt f = 0; f < Nfd1; f++) { 1313 // construct all formDegree+1 forms that start with dx_(dim - d) /\ ... 1314 form_atoms[0] = dim - d; 1315 PetscCall(PetscDTEnumSubset(d - 1, formDegree, f, &form_atoms[1])); 1316 for (PetscInt i = 0; i < formDegree; i++) form_atoms[1 + i] += form_atoms[0] + 1; 1317 PetscInt f_ind; // index of the resulting form 1318 PetscCall(PetscDTSubsetIndex(dim, formDegree + 1, form_atoms, &f_ind)); 1319 PetscReal *p_f = &p[total++ * Nf * Nk * npoints]; 1320 for (PetscInt nz = 0; nz < nnz; nz++) { 1321 PetscInt i = pattern[nz][0]; // formDegree component 1322 PetscInt j = pattern[nz][1]; // (formDegree + 1) component 1323 PetscInt v = pattern[nz][2]; // coordinate component 1324 PetscReal scale = v < 0 ? -1. : 1.; 1325 1326 i = formNegative ? (Nf - 1 - i) : i; 1327 scale = (formNegative && (i & 1)) ? -scale : scale; 1328 v = v < 0 ? -(v + 1) : v; 1329 if (j != f_ind) continue; 1330 PetscReal *p_i = &p_f[i * Nk * npoints]; 1331 for (PetscInt jet = 0; jet < Nk; jet++) { 1332 const PetscReal *h_jet = &h_s[jet * npoints]; 1333 PetscReal *p_jet = &p_i[jet * npoints]; 1334 1335 for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * h_jet[pt] * (points[pt * dim + v] - centroid); 1336 PetscCall(PetscDTIndexToGradedOrder(dim, jet, deriv)); 1337 deriv[v]++; 1338 PetscReal mult = deriv[v]; 1339 PetscInt l; 1340 PetscCall(PetscDTGradedOrderToIndex(dim, deriv, &l)); 1341 if (l >= Nk) continue; 1342 p_jet = &p_i[l * npoints]; 1343 for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * mult * h_jet[pt]; 1344 deriv[v]--; 1345 } 1346 } 1347 } 1348 } 1349 } 1350 PetscCheck(total == Nbpt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrectly counted P trimmed polynomials"); 1351 PetscCall(PetscFree(deriv)); 1352 PetscCall(PetscFree(pattern)); 1353 PetscCall(PetscFree(form_atoms)); 1354 PetscCall(PetscFree(p_scalar)); 1355 PetscFunctionReturn(PETSC_SUCCESS); 1356 } 1357 1358 /*@ 1359 PetscDTPTrimmedEvalJet - Evaluate the jet (function and derivatives) of a basis of the trimmed polynomial k-forms up to 1360 a given degree. 1361 1362 Input Parameters: 1363 + dim - the number of variables in the multivariate polynomials 1364 . npoints - the number of points to evaluate the polynomials at 1365 . points - [npoints x dim] array of point coordinates 1366 . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space to evaluate. 1367 There are ((dim + degree) choose (dim + formDegree)) x ((degree + formDegree - 1) choose (formDegree)) polynomials in this space. 1368 (You can use `PetscDTPTrimmedSize()` to compute this size.) 1369 . formDegree - the degree of the form 1370 - jetDegree - the maximum order partial derivative to evaluate in the jet. There are ((dim + jetDegree) choose dim) partial derivatives 1371 in the jet. Choosing jetDegree = 0 means to evaluate just the function and no derivatives 1372 1373 Output Parameter: 1374 . p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is 1375 `PetscDTPTrimmedSize()` x ((dim + formDegree) choose dim) x ((dim + k) choose dim) x npoints, 1376 which also describes the order of the dimensions of this 1377 four-dimensional array: 1378 the first (slowest varying) dimension is basis function index; 1379 the second dimension is component of the form; 1380 the third dimension is jet index; 1381 the fourth (fastest varying) dimension is the index of the evaluation point. 1382 1383 Level: advanced 1384 1385 Notes: 1386 The ordering of the basis functions is not graded, so the basis functions are not nested by degree like `PetscDTPKDEvalJet()`. 1387 The basis functions are not an L2-orthonormal basis on any particular domain. 1388 1389 The implementation is based on the description of the trimmed polynomials up to degree r as 1390 the direct sum of polynomials up to degree (r-1) and the Koszul differential applied to 1391 homogeneous polynomials of degree (r-1). 1392 1393 .seealso: `PetscDTPKDEvalJet()`, `PetscDTPTrimmedSize()` 1394 @*/ 1395 PetscErrorCode PetscDTPTrimmedEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[]) 1396 { 1397 PetscFunctionBegin; 1398 PetscCall(PetscDTPTrimmedEvalJet_Internal(dim, npoints, points, degree, formDegree, jetDegree, p)); 1399 PetscFunctionReturn(PETSC_SUCCESS); 1400 } 1401 1402 /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V 1403 * with lds n; diag and subdiag are overwritten */ 1404 static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[], PetscReal eigs[], PetscScalar V[]) 1405 { 1406 char jobz = 'V'; /* eigenvalues and eigenvectors */ 1407 char range = 'A'; /* all eigenvalues will be found */ 1408 PetscReal VL = 0.; /* ignored because range is 'A' */ 1409 PetscReal VU = 0.; /* ignored because range is 'A' */ 1410 PetscBLASInt IL = 0; /* ignored because range is 'A' */ 1411 PetscBLASInt IU = 0; /* ignored because range is 'A' */ 1412 PetscReal abstol = 0.; /* unused */ 1413 PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */ 1414 PetscBLASInt *isuppz; 1415 PetscBLASInt lwork, liwork; 1416 PetscReal workquery; 1417 PetscBLASInt iworkquery; 1418 PetscBLASInt *iwork; 1419 PetscBLASInt info; 1420 PetscReal *work = NULL; 1421 1422 PetscFunctionBegin; 1423 #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG) 1424 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found"); 1425 #endif 1426 PetscCall(PetscBLASIntCast(n, &bn)); 1427 PetscCall(PetscBLASIntCast(n, &ldz)); 1428 #if !defined(PETSC_MISSING_LAPACK_STEGR) 1429 PetscCall(PetscMalloc1(2 * n, &isuppz)); 1430 lwork = -1; 1431 liwork = -1; 1432 PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, &workquery, &lwork, &iworkquery, &liwork, &info)); 1433 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEGR error"); 1434 lwork = (PetscBLASInt)workquery; 1435 liwork = (PetscBLASInt)iworkquery; 1436 PetscCall(PetscMalloc2(lwork, &work, liwork, &iwork)); 1437 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 1438 PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, work, &lwork, iwork, &liwork, &info)); 1439 PetscCall(PetscFPTrapPop()); 1440 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEGR error"); 1441 PetscCall(PetscFree2(work, iwork)); 1442 PetscCall(PetscFree(isuppz)); 1443 #elif !defined(PETSC_MISSING_LAPACK_STEQR) 1444 jobz = 'I'; /* Compute eigenvalues and eigenvectors of the 1445 tridiagonal matrix. Z is initialized to the identity 1446 matrix. */ 1447 PetscCall(PetscMalloc1(PetscMax(1, 2 * n - 2), &work)); 1448 PetscCallBLAS("LAPACKsteqr", LAPACKsteqr_("I", &bn, diag, subdiag, V, &ldz, work, &info)); 1449 PetscCall(PetscFPTrapPop()); 1450 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEQR error"); 1451 PetscCall(PetscFree(work)); 1452 PetscCall(PetscArraycpy(eigs, diag, n)); 1453 #endif 1454 PetscFunctionReturn(PETSC_SUCCESS); 1455 } 1456 1457 /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi 1458 * quadrature rules on the interval [-1, 1] */ 1459 static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw) 1460 { 1461 PetscReal twoab1; 1462 PetscInt m = n - 2; 1463 PetscReal a = alpha + 1.; 1464 PetscReal b = beta + 1.; 1465 PetscReal gra, grb; 1466 1467 PetscFunctionBegin; 1468 twoab1 = PetscPowReal(2., a + b - 1.); 1469 #if defined(PETSC_HAVE_LGAMMA) 1470 grb = PetscExpReal(2. * PetscLGamma(b + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + a + 1.) - (PetscLGamma(m + b + 1) + PetscLGamma(m + a + b + 1.))); 1471 gra = PetscExpReal(2. * PetscLGamma(a + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + b + 1.) - (PetscLGamma(m + a + 1) + PetscLGamma(m + a + b + 1.))); 1472 #else 1473 { 1474 PetscInt alphai = (PetscInt)alpha; 1475 PetscInt betai = (PetscInt)beta; 1476 1477 if ((PetscReal)alphai == alpha && (PetscReal)betai == beta) { 1478 PetscReal binom1, binom2; 1479 1480 PetscCall(PetscDTBinomial(m + b, b, &binom1)); 1481 PetscCall(PetscDTBinomial(m + a + b, b, &binom2)); 1482 grb = 1. / (binom1 * binom2); 1483 PetscCall(PetscDTBinomial(m + a, a, &binom1)); 1484 PetscCall(PetscDTBinomial(m + a + b, a, &binom2)); 1485 gra = 1. / (binom1 * binom2); 1486 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable."); 1487 } 1488 #endif 1489 *leftw = twoab1 * grb / b; 1490 *rightw = twoab1 * gra / a; 1491 PetscFunctionReturn(PETSC_SUCCESS); 1492 } 1493 1494 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x. 1495 Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */ 1496 static inline PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 1497 { 1498 PetscReal pn1, pn2; 1499 PetscReal cnm1, cnm1x, cnm2; 1500 PetscInt k; 1501 1502 PetscFunctionBegin; 1503 if (!n) { 1504 *P = 1.0; 1505 PetscFunctionReturn(PETSC_SUCCESS); 1506 } 1507 PetscDTJacobiRecurrence_Internal(1, a, b, cnm1, cnm1x, cnm2); 1508 pn2 = 1.; 1509 pn1 = cnm1 + cnm1x * x; 1510 if (n == 1) { 1511 *P = pn1; 1512 PetscFunctionReturn(PETSC_SUCCESS); 1513 } 1514 *P = 0.0; 1515 for (k = 2; k < n + 1; ++k) { 1516 PetscDTJacobiRecurrence_Internal(k, a, b, cnm1, cnm1x, cnm2); 1517 1518 *P = (cnm1 + cnm1x * x) * pn1 - cnm2 * pn2; 1519 pn2 = pn1; 1520 pn1 = *P; 1521 } 1522 PetscFunctionReturn(PETSC_SUCCESS); 1523 } 1524 1525 /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */ 1526 static inline PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P) 1527 { 1528 PetscReal nP; 1529 PetscInt i; 1530 1531 PetscFunctionBegin; 1532 *P = 0.0; 1533 if (k > n) PetscFunctionReturn(PETSC_SUCCESS); 1534 PetscCall(PetscDTComputeJacobi(a + k, b + k, n - k, x, &nP)); 1535 for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5; 1536 *P = nP; 1537 PetscFunctionReturn(PETSC_SUCCESS); 1538 } 1539 1540 static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[]) 1541 { 1542 PetscInt maxIter = 100; 1543 PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON)); 1544 PetscReal a1, a6, gf; 1545 PetscInt k; 1546 1547 PetscFunctionBegin; 1548 1549 a1 = PetscPowReal(2.0, a + b + 1); 1550 #if defined(PETSC_HAVE_LGAMMA) 1551 { 1552 PetscReal a2, a3, a4, a5; 1553 a2 = PetscLGamma(a + npoints + 1); 1554 a3 = PetscLGamma(b + npoints + 1); 1555 a4 = PetscLGamma(a + b + npoints + 1); 1556 a5 = PetscLGamma(npoints + 1); 1557 gf = PetscExpReal(a2 + a3 - (a4 + a5)); 1558 } 1559 #else 1560 { 1561 PetscInt ia, ib; 1562 1563 ia = (PetscInt)a; 1564 ib = (PetscInt)b; 1565 gf = 1.; 1566 if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */ 1567 for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k); 1568 } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */ 1569 for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k); 1570 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable."); 1571 } 1572 #endif 1573 1574 a6 = a1 * gf; 1575 /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses. 1576 Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */ 1577 for (k = 0; k < npoints; ++k) { 1578 PetscReal r = PetscCosReal(PETSC_PI * (1. - (4. * k + 3. + 2. * b) / (4. * npoints + 2. * (a + b + 1.)))), dP; 1579 PetscInt j; 1580 1581 if (k > 0) r = 0.5 * (r + x[k - 1]); 1582 for (j = 0; j < maxIter; ++j) { 1583 PetscReal s = 0.0, delta, f, fp; 1584 PetscInt i; 1585 1586 for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]); 1587 PetscCall(PetscDTComputeJacobi(a, b, npoints, r, &f)); 1588 PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp)); 1589 delta = f / (fp - f * s); 1590 r = r - delta; 1591 if (PetscAbsReal(delta) < eps) break; 1592 } 1593 x[k] = r; 1594 PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP)); 1595 w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP); 1596 } 1597 PetscFunctionReturn(PETSC_SUCCESS); 1598 } 1599 1600 /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi 1601 * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */ 1602 static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s) 1603 { 1604 PetscInt i; 1605 1606 PetscFunctionBegin; 1607 for (i = 0; i < nPoints; i++) { 1608 PetscReal A, B, C; 1609 1610 PetscDTJacobiRecurrence_Internal(i + 1, a, b, A, B, C); 1611 d[i] = -A / B; 1612 if (i) s[i - 1] *= C / B; 1613 if (i < nPoints - 1) s[i] = 1. / B; 1614 } 1615 PetscFunctionReturn(PETSC_SUCCESS); 1616 } 1617 1618 static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w) 1619 { 1620 PetscReal mu0; 1621 PetscReal ga, gb, gab; 1622 PetscInt i; 1623 1624 PetscFunctionBegin; 1625 PetscCall(PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite)); 1626 1627 #if defined(PETSC_HAVE_TGAMMA) 1628 ga = PetscTGamma(a + 1); 1629 gb = PetscTGamma(b + 1); 1630 gab = PetscTGamma(a + b + 2); 1631 #else 1632 { 1633 PetscInt ia, ib; 1634 1635 ia = (PetscInt)a; 1636 ib = (PetscInt)b; 1637 if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */ 1638 PetscCall(PetscDTFactorial(ia, &ga)); 1639 PetscCall(PetscDTFactorial(ib, &gb)); 1640 PetscCall(PetscDTFactorial(ia + ib + 1, &gb)); 1641 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "tgamma() - math routine is unavailable."); 1642 } 1643 #endif 1644 mu0 = PetscPowReal(2., a + b + 1.) * ga * gb / gab; 1645 1646 #if defined(PETSCDTGAUSSIANQUADRATURE_EIG) 1647 { 1648 PetscReal *diag, *subdiag; 1649 PetscScalar *V; 1650 1651 PetscCall(PetscMalloc2(npoints, &diag, npoints, &subdiag)); 1652 PetscCall(PetscMalloc1(npoints * npoints, &V)); 1653 PetscCall(PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag)); 1654 for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]); 1655 PetscCall(PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V)); 1656 for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0; 1657 PetscCall(PetscFree(V)); 1658 PetscCall(PetscFree2(diag, subdiag)); 1659 } 1660 #else 1661 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found"); 1662 #endif 1663 { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the 1664 eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that 1665 the eigenvalues are sorted */ 1666 PetscBool sorted; 1667 1668 PetscCall(PetscSortedReal(npoints, x, &sorted)); 1669 if (!sorted) { 1670 PetscInt *order, i; 1671 PetscReal *tmp; 1672 1673 PetscCall(PetscMalloc2(npoints, &order, npoints, &tmp)); 1674 for (i = 0; i < npoints; i++) order[i] = i; 1675 PetscCall(PetscSortRealWithPermutation(npoints, x, order)); 1676 PetscCall(PetscArraycpy(tmp, x, npoints)); 1677 for (i = 0; i < npoints; i++) x[i] = tmp[order[i]]; 1678 PetscCall(PetscArraycpy(tmp, w, npoints)); 1679 for (i = 0; i < npoints; i++) w[i] = tmp[order[i]]; 1680 PetscCall(PetscFree2(order, tmp)); 1681 } 1682 } 1683 PetscFunctionReturn(PETSC_SUCCESS); 1684 } 1685 1686 static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1687 { 1688 PetscFunctionBegin; 1689 PetscCheck(npoints >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive"); 1690 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1691 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1."); 1692 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1."); 1693 1694 if (newton) PetscCall(PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w)); 1695 else PetscCall(PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w)); 1696 if (alpha == beta) { /* symmetrize */ 1697 PetscInt i; 1698 for (i = 0; i < (npoints + 1) / 2; i++) { 1699 PetscInt j = npoints - 1 - i; 1700 PetscReal xi = x[i]; 1701 PetscReal xj = x[j]; 1702 PetscReal wi = w[i]; 1703 PetscReal wj = w[j]; 1704 1705 x[i] = (xi - xj) / 2.; 1706 x[j] = (xj - xi) / 2.; 1707 w[i] = w[j] = (wi + wj) / 2.; 1708 } 1709 } 1710 PetscFunctionReturn(PETSC_SUCCESS); 1711 } 1712 1713 /*@ 1714 PetscDTGaussJacobiQuadrature - quadrature for the interval [a, b] with the weight function 1715 $(x-a)^\alpha (x-b)^\beta$. 1716 1717 Not Collective 1718 1719 Input Parameters: 1720 + npoints - the number of points in the quadrature rule 1721 . a - the left endpoint of the interval 1722 . b - the right endpoint of the interval 1723 . alpha - the left exponent 1724 - beta - the right exponent 1725 1726 Output Parameters: 1727 + x - array of length `npoints`, the locations of the quadrature points 1728 - w - array of length `npoints`, the weights of the quadrature points 1729 1730 Level: intermediate 1731 1732 Note: 1733 This quadrature rule is exact for polynomials up to degree 2*npoints - 1. 1734 1735 .seealso: `PetscDTGaussQuadrature()` 1736 @*/ 1737 PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1738 { 1739 PetscInt i; 1740 1741 PetscFunctionBegin; 1742 PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal)); 1743 if (a != -1. || b != 1.) { /* shift */ 1744 for (i = 0; i < npoints; i++) { 1745 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1746 w[i] *= (b - a) / 2.; 1747 } 1748 } 1749 PetscFunctionReturn(PETSC_SUCCESS); 1750 } 1751 1752 static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1753 { 1754 PetscInt i; 1755 1756 PetscFunctionBegin; 1757 PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive"); 1758 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1759 PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1."); 1760 PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1."); 1761 1762 x[0] = -1.; 1763 x[npoints - 1] = 1.; 1764 if (npoints > 2) PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints - 2, alpha + 1., beta + 1., &x[1], &w[1], newton)); 1765 for (i = 1; i < npoints - 1; i++) w[i] /= (1. - x[i] * x[i]); 1766 PetscCall(PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints - 1])); 1767 PetscFunctionReturn(PETSC_SUCCESS); 1768 } 1769 1770 /*@ 1771 PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval [a, b] with the weight function 1772 $(x-a)^\alpha (x-b)^\beta$, with endpoints a and b included as quadrature points. 1773 1774 Not Collective 1775 1776 Input Parameters: 1777 + npoints - the number of points in the quadrature rule 1778 . a - the left endpoint of the interval 1779 . b - the right endpoint of the interval 1780 . alpha - the left exponent 1781 - beta - the right exponent 1782 1783 Output Parameters: 1784 + x - array of length `npoints`, the locations of the quadrature points 1785 - w - array of length `npoints`, the weights of the quadrature points 1786 1787 Level: intermediate 1788 1789 Note: 1790 This quadrature rule is exact for polynomials up to degree 2*npoints - 3. 1791 1792 .seealso: `PetscDTGaussJacobiQuadrature()` 1793 @*/ 1794 PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1795 { 1796 PetscInt i; 1797 1798 PetscFunctionBegin; 1799 PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal)); 1800 if (a != -1. || b != 1.) { /* shift */ 1801 for (i = 0; i < npoints; i++) { 1802 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1803 w[i] *= (b - a) / 2.; 1804 } 1805 } 1806 PetscFunctionReturn(PETSC_SUCCESS); 1807 } 1808 1809 /*@ 1810 PetscDTGaussQuadrature - create Gauss-Legendre quadrature 1811 1812 Not Collective 1813 1814 Input Parameters: 1815 + npoints - number of points 1816 . a - left end of interval (often-1) 1817 - b - right end of interval (often +1) 1818 1819 Output Parameters: 1820 + x - quadrature points 1821 - w - quadrature weights 1822 1823 Level: intermediate 1824 1825 References: 1826 . * - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969. 1827 1828 .seealso: `PetscDTLegendreEval()`, `PetscDTGaussJacobiQuadrature()` 1829 @*/ 1830 PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w) 1831 { 1832 PetscInt i; 1833 1834 PetscFunctionBegin; 1835 PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal)); 1836 if (a != -1. || b != 1.) { /* shift */ 1837 for (i = 0; i < npoints; i++) { 1838 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1839 w[i] *= (b - a) / 2.; 1840 } 1841 } 1842 PetscFunctionReturn(PETSC_SUCCESS); 1843 } 1844 1845 /*@C 1846 PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre 1847 nodes of a given size on the domain [-1,1] 1848 1849 Not Collective 1850 1851 Input Parameters: 1852 + n - number of grid nodes 1853 - type - `PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA` or `PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON` 1854 1855 Output Parameters: 1856 + x - quadrature points 1857 - w - quadrature weights 1858 1859 Level: intermediate 1860 1861 Notes: 1862 For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not 1863 close enough to the desired solution 1864 1865 These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes 1866 1867 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 1868 1869 .seealso: `PetscDTGaussQuadrature()`, `PetscGaussLobattoLegendreCreateType` 1870 1871 @*/ 1872 PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints, PetscGaussLobattoLegendreCreateType type, PetscReal *x, PetscReal *w) 1873 { 1874 PetscBool newton; 1875 1876 PetscFunctionBegin; 1877 PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide at least 2 grid points per element"); 1878 newton = (PetscBool)(type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON); 1879 PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton)); 1880 PetscFunctionReturn(PETSC_SUCCESS); 1881 } 1882 1883 /*@ 1884 PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature 1885 1886 Not Collective 1887 1888 Input Parameters: 1889 + dim - The spatial dimension 1890 . Nc - The number of components 1891 . npoints - number of points in one dimension 1892 . a - left end of interval (often-1) 1893 - b - right end of interval (often +1) 1894 1895 Output Parameter: 1896 . q - A `PetscQuadrature` object 1897 1898 Level: intermediate 1899 1900 .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()` 1901 @*/ 1902 PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 1903 { 1904 DMPolytopeType ct; 1905 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints * PetscSqr(npoints) : PetscSqr(npoints) : npoints; 1906 PetscReal *x, *w, *xw, *ww; 1907 1908 PetscFunctionBegin; 1909 PetscCall(PetscMalloc1(totpoints * dim, &x)); 1910 PetscCall(PetscMalloc1(totpoints * Nc, &w)); 1911 /* Set up the Golub-Welsch system */ 1912 switch (dim) { 1913 case 0: 1914 ct = DM_POLYTOPE_POINT; 1915 PetscCall(PetscFree(x)); 1916 PetscCall(PetscFree(w)); 1917 PetscCall(PetscMalloc1(1, &x)); 1918 PetscCall(PetscMalloc1(Nc, &w)); 1919 x[0] = 0.0; 1920 for (PetscInt c = 0; c < Nc; ++c) w[c] = 1.0; 1921 break; 1922 case 1: 1923 ct = DM_POLYTOPE_SEGMENT; 1924 PetscCall(PetscMalloc1(npoints, &ww)); 1925 PetscCall(PetscDTGaussQuadrature(npoints, a, b, x, ww)); 1926 for (PetscInt i = 0; i < npoints; ++i) 1927 for (PetscInt c = 0; c < Nc; ++c) w[i * Nc + c] = ww[i]; 1928 PetscCall(PetscFree(ww)); 1929 break; 1930 case 2: 1931 ct = DM_POLYTOPE_QUADRILATERAL; 1932 PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww)); 1933 PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww)); 1934 for (PetscInt i = 0; i < npoints; ++i) { 1935 for (PetscInt j = 0; j < npoints; ++j) { 1936 x[(i * npoints + j) * dim + 0] = xw[i]; 1937 x[(i * npoints + j) * dim + 1] = xw[j]; 1938 for (PetscInt c = 0; c < Nc; ++c) w[(i * npoints + j) * Nc + c] = ww[i] * ww[j]; 1939 } 1940 } 1941 PetscCall(PetscFree2(xw, ww)); 1942 break; 1943 case 3: 1944 ct = DM_POLYTOPE_HEXAHEDRON; 1945 PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww)); 1946 PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww)); 1947 for (PetscInt i = 0; i < npoints; ++i) { 1948 for (PetscInt j = 0; j < npoints; ++j) { 1949 for (PetscInt k = 0; k < npoints; ++k) { 1950 x[((i * npoints + j) * npoints + k) * dim + 0] = xw[i]; 1951 x[((i * npoints + j) * npoints + k) * dim + 1] = xw[j]; 1952 x[((i * npoints + j) * npoints + k) * dim + 2] = xw[k]; 1953 for (PetscInt c = 0; c < Nc; ++c) w[((i * npoints + j) * npoints + k) * Nc + c] = ww[i] * ww[j] * ww[k]; 1954 } 1955 } 1956 } 1957 PetscCall(PetscFree2(xw, ww)); 1958 break; 1959 default: 1960 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %" PetscInt_FMT, dim); 1961 } 1962 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 1963 PetscCall(PetscQuadratureSetCellType(*q, ct)); 1964 PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1)); 1965 PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w)); 1966 PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "GaussTensor")); 1967 PetscFunctionReturn(PETSC_SUCCESS); 1968 } 1969 1970 /*@ 1971 PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex 1972 1973 Not Collective 1974 1975 Input Parameters: 1976 + dim - The simplex dimension 1977 . Nc - The number of components 1978 . npoints - The number of points in one dimension 1979 . a - left end of interval (often-1) 1980 - b - right end of interval (often +1) 1981 1982 Output Parameter: 1983 . q - A `PetscQuadrature` object 1984 1985 Level: intermediate 1986 1987 Note: 1988 For `dim` == 1, this is Gauss-Legendre quadrature 1989 1990 References: 1991 . * - Karniadakis and Sherwin. FIAT 1992 1993 .seealso: `PetscDTGaussTensorQuadrature()`, `PetscDTGaussQuadrature()` 1994 @*/ 1995 PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 1996 { 1997 DMPolytopeType ct; 1998 PetscInt totpoints; 1999 PetscReal *p1, *w1; 2000 PetscReal *x, *w; 2001 2002 PetscFunctionBegin; 2003 PetscCheck(!(a != -1.0) && !(b != 1.0), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now"); 2004 switch (dim) { 2005 case 0: 2006 ct = DM_POLYTOPE_POINT; 2007 break; 2008 case 1: 2009 ct = DM_POLYTOPE_SEGMENT; 2010 break; 2011 case 2: 2012 ct = DM_POLYTOPE_TRIANGLE; 2013 break; 2014 case 3: 2015 ct = DM_POLYTOPE_TETRAHEDRON; 2016 break; 2017 default: 2018 ct = DM_POLYTOPE_UNKNOWN; 2019 } 2020 totpoints = 1; 2021 for (PetscInt i = 0; i < dim; ++i) totpoints *= npoints; 2022 PetscCall(PetscMalloc1(totpoints * dim, &x)); 2023 PetscCall(PetscMalloc1(totpoints * Nc, &w)); 2024 PetscCall(PetscMalloc2(npoints, &p1, npoints, &w1)); 2025 for (PetscInt i = 0; i < totpoints * Nc; ++i) w[i] = 1.; 2026 for (PetscInt i = 0, totprev = 1, totrem = totpoints / npoints; i < dim; ++i) { 2027 PetscReal mul; 2028 2029 mul = PetscPowReal(2., -i); 2030 PetscCall(PetscDTGaussJacobiQuadrature(npoints, -1., 1., i, 0.0, p1, w1)); 2031 for (PetscInt pt = 0, l = 0; l < totprev; l++) { 2032 for (PetscInt j = 0; j < npoints; j++) { 2033 for (PetscInt m = 0; m < totrem; m++, pt++) { 2034 for (PetscInt k = 0; k < i; k++) x[pt * dim + k] = (x[pt * dim + k] + 1.) * (1. - p1[j]) * 0.5 - 1.; 2035 x[pt * dim + i] = p1[j]; 2036 for (PetscInt c = 0; c < Nc; c++) w[pt * Nc + c] *= mul * w1[j]; 2037 } 2038 } 2039 } 2040 totprev *= npoints; 2041 totrem /= npoints; 2042 } 2043 PetscCall(PetscFree2(p1, w1)); 2044 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 2045 PetscCall(PetscQuadratureSetCellType(*q, ct)); 2046 PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1)); 2047 PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w)); 2048 PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "StroudConical")); 2049 PetscFunctionReturn(PETSC_SUCCESS); 2050 } 2051 2052 static PetscBool MinSymTriQuadCite = PETSC_FALSE; 2053 const char MinSymTriQuadCitation[] = "@article{WitherdenVincent2015,\n" 2054 " title = {On the identification of symmetric quadrature rules for finite element methods},\n" 2055 " journal = {Computers & Mathematics with Applications},\n" 2056 " volume = {69},\n" 2057 " number = {10},\n" 2058 " pages = {1232-1241},\n" 2059 " year = {2015},\n" 2060 " issn = {0898-1221},\n" 2061 " doi = {10.1016/j.camwa.2015.03.017},\n" 2062 " url = {https://www.sciencedirect.com/science/article/pii/S0898122115001224},\n" 2063 " author = {F.D. Witherden and P.E. Vincent},\n" 2064 "}\n"; 2065 2066 #include "petscdttriquadrules.h" 2067 2068 static PetscBool MinSymTetQuadCite = PETSC_FALSE; 2069 const char MinSymTetQuadCitation[] = "@article{JaskowiecSukumar2021\n" 2070 " author = {Jaskowiec, Jan and Sukumar, N.},\n" 2071 " title = {High-order symmetric cubature rules for tetrahedra and pyramids},\n" 2072 " journal = {International Journal for Numerical Methods in Engineering},\n" 2073 " volume = {122},\n" 2074 " number = {1},\n" 2075 " pages = {148-171},\n" 2076 " doi = {10.1002/nme.6528},\n" 2077 " url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6528},\n" 2078 " eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6528},\n" 2079 " year = {2021}\n" 2080 "}\n"; 2081 2082 #include "petscdttetquadrules.h" 2083 2084 // https://en.wikipedia.org/wiki/Partition_(number_theory) 2085 static PetscErrorCode PetscDTPartitionNumber(PetscInt n, PetscInt *p) 2086 { 2087 // sequence A000041 in the OEIS 2088 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}; 2089 PetscInt tabulated_max = PETSC_STATIC_ARRAY_LENGTH(partition) - 1; 2090 2091 PetscFunctionBegin; 2092 PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Partition number not defined for negative number %" PetscInt_FMT, n); 2093 // not implementing the pentagonal number recurrence, we don't need partition numbers for n that high 2094 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); 2095 *p = partition[n]; 2096 PetscFunctionReturn(PETSC_SUCCESS); 2097 } 2098 2099 /*@ 2100 PetscDTSimplexQuadrature - Create a quadrature rule for a simplex that exactly integrates polynomials up to a given degree. 2101 2102 Not Collective 2103 2104 Input Parameters: 2105 + dim - The spatial dimension of the simplex (1 = segment, 2 = triangle, 3 = tetrahedron) 2106 . degree - The largest polynomial degree that is required to be integrated exactly 2107 - type - left end of interval (often-1) 2108 2109 Output Parameter: 2110 . quad - A `PetscQuadrature` object for integration over the biunit simplex 2111 (defined by the bounds $x_i >= -1$ and $\sum_i x_i <= 2 - d$) that is exact for 2112 polynomials up to the given degree 2113 2114 Level: intermediate 2115 2116 .seealso: `PetscDTSimplexQuadratureType`, `PetscDTGaussQuadrature()`, `PetscDTStroudCononicalQuadrature()`, `PetscQuadrature` 2117 @*/ 2118 PetscErrorCode PetscDTSimplexQuadrature(PetscInt dim, PetscInt degree, PetscDTSimplexQuadratureType type, PetscQuadrature *quad) 2119 { 2120 PetscDTSimplexQuadratureType orig_type = type; 2121 2122 PetscFunctionBegin; 2123 PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative dimension %" PetscInt_FMT, dim); 2124 PetscCheck(degree >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %" PetscInt_FMT, degree); 2125 if (type == PETSCDTSIMPLEXQUAD_DEFAULT) type = PETSCDTSIMPLEXQUAD_MINSYM; 2126 if (type == PETSCDTSIMPLEXQUAD_CONIC || dim < 2) { 2127 PetscInt points_per_dim = (degree + 2) / 2; // ceil((degree + 1) / 2); 2128 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, points_per_dim, -1, 1, quad)); 2129 } else { 2130 DMPolytopeType ct; 2131 PetscInt n = dim + 1; 2132 PetscInt fact = 1; 2133 PetscInt *part, *perm; 2134 PetscInt p = 0; 2135 PetscInt max_degree; 2136 const PetscInt *nodes_per_type = NULL; 2137 const PetscInt *all_num_full_nodes = NULL; 2138 const PetscReal **weights_list = NULL; 2139 const PetscReal **compact_nodes_list = NULL; 2140 const char *citation = NULL; 2141 PetscBool *cited = NULL; 2142 2143 switch (dim) { 2144 case 0: 2145 ct = DM_POLYTOPE_POINT; 2146 break; 2147 case 1: 2148 ct = DM_POLYTOPE_SEGMENT; 2149 break; 2150 case 2: 2151 ct = DM_POLYTOPE_TRIANGLE; 2152 break; 2153 case 3: 2154 ct = DM_POLYTOPE_TETRAHEDRON; 2155 break; 2156 default: 2157 ct = DM_POLYTOPE_UNKNOWN; 2158 } 2159 switch (dim) { 2160 case 2: 2161 cited = &MinSymTriQuadCite; 2162 citation = MinSymTriQuadCitation; 2163 max_degree = PetscDTWVTriQuad_max_degree; 2164 nodes_per_type = PetscDTWVTriQuad_num_orbits; 2165 all_num_full_nodes = PetscDTWVTriQuad_num_nodes; 2166 weights_list = PetscDTWVTriQuad_weights; 2167 compact_nodes_list = PetscDTWVTriQuad_orbits; 2168 break; 2169 case 3: 2170 cited = &MinSymTetQuadCite; 2171 citation = MinSymTetQuadCitation; 2172 max_degree = PetscDTJSTetQuad_max_degree; 2173 nodes_per_type = PetscDTJSTetQuad_num_orbits; 2174 all_num_full_nodes = PetscDTJSTetQuad_num_nodes; 2175 weights_list = PetscDTJSTetQuad_weights; 2176 compact_nodes_list = PetscDTJSTetQuad_orbits; 2177 break; 2178 default: 2179 max_degree = -1; 2180 break; 2181 } 2182 2183 if (degree > max_degree) { 2184 if (orig_type == PETSCDTSIMPLEXQUAD_DEFAULT) { 2185 // fall back to conic 2186 PetscCall(PetscDTSimplexQuadrature(dim, degree, PETSCDTSIMPLEXQUAD_CONIC, quad)); 2187 PetscFunctionReturn(PETSC_SUCCESS); 2188 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Minimal symmetric quadrature for dim %" PetscInt_FMT ", degree %" PetscInt_FMT " unsupported", dim, degree); 2189 } 2190 2191 PetscCall(PetscCitationsRegister(citation, cited)); 2192 2193 PetscCall(PetscDTPartitionNumber(n, &p)); 2194 for (PetscInt d = 2; d <= n; d++) fact *= d; 2195 2196 PetscInt num_full_nodes = all_num_full_nodes[degree]; 2197 const PetscReal *all_compact_nodes = compact_nodes_list[degree]; 2198 const PetscReal *all_compact_weights = weights_list[degree]; 2199 nodes_per_type = &nodes_per_type[p * degree]; 2200 2201 PetscReal *points; 2202 PetscReal *counts; 2203 PetscReal *weights; 2204 PetscReal *bary_to_biunit; // row-major transformation of barycentric coordinate to biunit 2205 PetscQuadrature q; 2206 2207 // compute the transformation 2208 PetscCall(PetscMalloc1(n * dim, &bary_to_biunit)); 2209 for (PetscInt d = 0; d < dim; d++) { 2210 for (PetscInt b = 0; b < n; b++) bary_to_biunit[d * n + b] = (d == b) ? 1.0 : -1.0; 2211 } 2212 2213 PetscCall(PetscMalloc3(n, &part, n, &perm, n, &counts)); 2214 PetscCall(PetscCalloc1(num_full_nodes * dim, &points)); 2215 PetscCall(PetscMalloc1(num_full_nodes, &weights)); 2216 2217 // (0, 0, ...) is the first partition lexicographically 2218 PetscCall(PetscArrayzero(part, n)); 2219 PetscCall(PetscArrayzero(counts, n)); 2220 counts[0] = n; 2221 2222 // for each partition 2223 for (PetscInt s = 0, node_offset = 0; s < p; s++) { 2224 PetscInt num_compact_coords = part[n - 1] + 1; 2225 2226 const PetscReal *compact_nodes = all_compact_nodes; 2227 const PetscReal *compact_weights = all_compact_weights; 2228 all_compact_nodes += num_compact_coords * nodes_per_type[s]; 2229 all_compact_weights += nodes_per_type[s]; 2230 2231 // for every permutation of the vertices 2232 for (PetscInt f = 0; f < fact; f++) { 2233 PetscCall(PetscDTEnumPerm(n, f, perm, NULL)); 2234 2235 // check if it is a valid permutation 2236 PetscInt digit; 2237 for (digit = 1; digit < n; digit++) { 2238 // skip permutations that would duplicate a node because it has a smaller symmetry group 2239 if (part[digit - 1] == part[digit] && perm[digit - 1] > perm[digit]) break; 2240 } 2241 if (digit < n) continue; 2242 2243 // create full nodes from this permutation of the compact nodes 2244 PetscReal *full_nodes = &points[node_offset * dim]; 2245 PetscReal *full_weights = &weights[node_offset]; 2246 2247 PetscCall(PetscArraycpy(full_weights, compact_weights, nodes_per_type[s])); 2248 for (PetscInt b = 0; b < n; b++) { 2249 for (PetscInt d = 0; d < dim; d++) { 2250 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]]; 2251 } 2252 } 2253 node_offset += nodes_per_type[s]; 2254 } 2255 2256 if (s < p - 1) { // Generate the next partition 2257 /* A partition is described by the number of coordinates that are in 2258 * each set of duplicates (counts) and redundantly by mapping each 2259 * index to its set of duplicates (part) 2260 * 2261 * Counts should always be in nonincreasing order 2262 * 2263 * We want to generate the partitions lexically by part, which means 2264 * finding the last index where count > 1 and reducing by 1. 2265 * 2266 * For the new counts beyond that index, we eagerly assign the remaining 2267 * capacity of the partition to smaller indices (ensures lexical ordering), 2268 * while respecting the nonincreasing invariant of the counts 2269 */ 2270 PetscInt last_digit = part[n - 1]; 2271 PetscInt last_digit_with_extra = last_digit; 2272 while (counts[last_digit_with_extra] == 1) last_digit_with_extra--; 2273 PetscInt limit = --counts[last_digit_with_extra]; 2274 PetscInt total_to_distribute = last_digit - last_digit_with_extra + 1; 2275 for (PetscInt digit = last_digit_with_extra + 1; digit < n; digit++) { 2276 counts[digit] = PetscMin(limit, total_to_distribute); 2277 total_to_distribute -= PetscMin(limit, total_to_distribute); 2278 } 2279 for (PetscInt digit = 0, offset = 0; digit < n; digit++) { 2280 PetscInt count = counts[digit]; 2281 for (PetscInt c = 0; c < count; c++) part[offset++] = digit; 2282 } 2283 } 2284 } 2285 PetscCall(PetscFree3(part, perm, counts)); 2286 PetscCall(PetscFree(bary_to_biunit)); 2287 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &q)); 2288 PetscCall(PetscQuadratureSetCellType(q, ct)); 2289 PetscCall(PetscQuadratureSetOrder(q, degree)); 2290 PetscCall(PetscQuadratureSetData(q, dim, 1, num_full_nodes, points, weights)); 2291 *quad = q; 2292 } 2293 PetscFunctionReturn(PETSC_SUCCESS); 2294 } 2295 2296 /*@ 2297 PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell 2298 2299 Not Collective 2300 2301 Input Parameters: 2302 + dim - The cell dimension 2303 . level - The number of points in one dimension, 2^l 2304 . a - left end of interval (often-1) 2305 - b - right end of interval (often +1) 2306 2307 Output Parameter: 2308 . q - A `PetscQuadrature` object 2309 2310 Level: intermediate 2311 2312 .seealso: `PetscDTGaussTensorQuadrature()`, `PetscQuadrature` 2313 @*/ 2314 PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q) 2315 { 2316 DMPolytopeType ct; 2317 const PetscInt p = 16; /* Digits of precision in the evaluation */ 2318 const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */ 2319 const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */ 2320 const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */ 2321 PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */ 2322 PetscReal wk = 0.5 * PETSC_PI; /* Quadrature weight at x_k */ 2323 PetscReal *x, *w; 2324 PetscInt K, k, npoints; 2325 2326 PetscFunctionBegin; 2327 PetscCheck(dim <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not yet implemented", dim); 2328 PetscCheck(level, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits"); 2329 switch (dim) { 2330 case 0: 2331 ct = DM_POLYTOPE_POINT; 2332 break; 2333 case 1: 2334 ct = DM_POLYTOPE_SEGMENT; 2335 break; 2336 case 2: 2337 ct = DM_POLYTOPE_QUADRILATERAL; 2338 break; 2339 case 3: 2340 ct = DM_POLYTOPE_HEXAHEDRON; 2341 break; 2342 default: 2343 ct = DM_POLYTOPE_UNKNOWN; 2344 } 2345 /* Find K such that the weights are < 32 digits of precision */ 2346 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))); 2347 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 2348 PetscCall(PetscQuadratureSetCellType(*q, ct)); 2349 PetscCall(PetscQuadratureSetOrder(*q, 2 * K + 1)); 2350 npoints = 2 * K - 1; 2351 PetscCall(PetscMalloc1(npoints * dim, &x)); 2352 PetscCall(PetscMalloc1(npoints, &w)); 2353 /* Center term */ 2354 x[0] = beta; 2355 w[0] = 0.5 * alpha * PETSC_PI; 2356 for (k = 1; k < K; ++k) { 2357 wk = 0.5 * alpha * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2358 xk = PetscTanhReal(0.5 * PETSC_PI * PetscSinhReal(k * h)); 2359 x[2 * k - 1] = -alpha * xk + beta; 2360 w[2 * k - 1] = wk; 2361 x[2 * k + 0] = alpha * xk + beta; 2362 w[2 * k + 0] = wk; 2363 } 2364 PetscCall(PetscQuadratureSetData(*q, dim, 1, npoints, x, w)); 2365 PetscFunctionReturn(PETSC_SUCCESS); 2366 } 2367 2368 PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2369 { 2370 const PetscInt p = 16; /* Digits of precision in the evaluation */ 2371 const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */ 2372 const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */ 2373 PetscReal h = 1.0; /* Step size, length between x_k */ 2374 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 2375 PetscReal osum = 0.0; /* Integral on last level */ 2376 PetscReal psum = 0.0; /* Integral on the level before the last level */ 2377 PetscReal sum; /* Integral on current level */ 2378 PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 2379 PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 2380 PetscReal wk; /* Quadrature weight at x_k */ 2381 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 2382 PetscInt d; /* Digits of precision in the integral */ 2383 2384 PetscFunctionBegin; 2385 PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 2386 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 2387 /* Center term */ 2388 func(&beta, ctx, &lval); 2389 sum = 0.5 * alpha * PETSC_PI * lval; 2390 /* */ 2391 do { 2392 PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4; 2393 PetscInt k = 1; 2394 2395 ++l; 2396 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */ 2397 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 2398 psum = osum; 2399 osum = sum; 2400 h *= 0.5; 2401 sum *= 0.5; 2402 do { 2403 wk = 0.5 * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2404 yk = 1.0 / (PetscExpReal(0.5 * PETSC_PI * PetscSinhReal(k * h)) * PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h))); 2405 lx = -alpha * (1.0 - yk) + beta; 2406 rx = alpha * (1.0 - yk) + beta; 2407 func(&lx, ctx, &lval); 2408 func(&rx, ctx, &rval); 2409 lterm = alpha * wk * lval; 2410 maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm); 2411 sum += lterm; 2412 rterm = alpha * wk * rval; 2413 maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm); 2414 sum += rterm; 2415 ++k; 2416 /* Only need to evaluate every other point on refined levels */ 2417 if (l != 1) ++k; 2418 } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */ 2419 2420 d1 = PetscLog10Real(PetscAbsReal(sum - osum)); 2421 d2 = PetscLog10Real(PetscAbsReal(sum - psum)); 2422 d3 = PetscLog10Real(maxTerm) - p; 2423 if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0; 2424 else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm))); 2425 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4))); 2426 } while (d < digits && l < 12); 2427 *sol = sum; 2428 PetscCall(PetscFPTrapPop()); 2429 PetscFunctionReturn(PETSC_SUCCESS); 2430 } 2431 2432 #if defined(PETSC_HAVE_MPFR) 2433 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2434 { 2435 const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */ 2436 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 2437 mpfr_t alpha; /* Half-width of the integration interval */ 2438 mpfr_t beta; /* Center of the integration interval */ 2439 mpfr_t h; /* Step size, length between x_k */ 2440 mpfr_t osum; /* Integral on last level */ 2441 mpfr_t psum; /* Integral on the level before the last level */ 2442 mpfr_t sum; /* Integral on current level */ 2443 mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 2444 mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 2445 mpfr_t wk; /* Quadrature weight at x_k */ 2446 PetscReal lval, rval, rtmp; /* Terms in the quadature sum to the left and right of 0 */ 2447 PetscInt d; /* Digits of precision in the integral */ 2448 mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp; 2449 2450 PetscFunctionBegin; 2451 PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 2452 /* Create high precision storage */ 2453 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); 2454 /* Initialization */ 2455 mpfr_set_d(alpha, 0.5 * (b - a), MPFR_RNDN); 2456 mpfr_set_d(beta, 0.5 * (b + a), MPFR_RNDN); 2457 mpfr_set_d(osum, 0.0, MPFR_RNDN); 2458 mpfr_set_d(psum, 0.0, MPFR_RNDN); 2459 mpfr_set_d(h, 1.0, MPFR_RNDN); 2460 mpfr_const_pi(pi2, MPFR_RNDN); 2461 mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN); 2462 /* Center term */ 2463 rtmp = 0.5 * (b + a); 2464 func(&rtmp, ctx, &lval); 2465 mpfr_set(sum, pi2, MPFR_RNDN); 2466 mpfr_mul(sum, sum, alpha, MPFR_RNDN); 2467 mpfr_mul_d(sum, sum, lval, MPFR_RNDN); 2468 /* */ 2469 do { 2470 PetscReal d1, d2, d3, d4; 2471 PetscInt k = 1; 2472 2473 ++l; 2474 mpfr_set_d(maxTerm, 0.0, MPFR_RNDN); 2475 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */ 2476 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 2477 mpfr_set(psum, osum, MPFR_RNDN); 2478 mpfr_set(osum, sum, MPFR_RNDN); 2479 mpfr_mul_d(h, h, 0.5, MPFR_RNDN); 2480 mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN); 2481 do { 2482 mpfr_set_si(kh, k, MPFR_RNDN); 2483 mpfr_mul(kh, kh, h, MPFR_RNDN); 2484 /* Weight */ 2485 mpfr_set(wk, h, MPFR_RNDN); 2486 mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN); 2487 mpfr_mul(msinh, msinh, pi2, MPFR_RNDN); 2488 mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN); 2489 mpfr_cosh(tmp, msinh, MPFR_RNDN); 2490 mpfr_sqr(tmp, tmp, MPFR_RNDN); 2491 mpfr_mul(wk, wk, mcosh, MPFR_RNDN); 2492 mpfr_div(wk, wk, tmp, MPFR_RNDN); 2493 /* Abscissa */ 2494 mpfr_set_d(yk, 1.0, MPFR_RNDZ); 2495 mpfr_cosh(tmp, msinh, MPFR_RNDN); 2496 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 2497 mpfr_exp(tmp, msinh, MPFR_RNDN); 2498 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 2499 /* Quadrature points */ 2500 mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ); 2501 mpfr_mul(lx, lx, alpha, MPFR_RNDU); 2502 mpfr_add(lx, lx, beta, MPFR_RNDU); 2503 mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ); 2504 mpfr_mul(rx, rx, alpha, MPFR_RNDD); 2505 mpfr_add(rx, rx, beta, MPFR_RNDD); 2506 /* Evaluation */ 2507 rtmp = mpfr_get_d(lx, MPFR_RNDU); 2508 func(&rtmp, ctx, &lval); 2509 rtmp = mpfr_get_d(rx, MPFR_RNDD); 2510 func(&rtmp, ctx, &rval); 2511 /* Update */ 2512 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 2513 mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN); 2514 mpfr_add(sum, sum, tmp, MPFR_RNDN); 2515 mpfr_abs(tmp, tmp, MPFR_RNDN); 2516 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 2517 mpfr_set(curTerm, tmp, MPFR_RNDN); 2518 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 2519 mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN); 2520 mpfr_add(sum, sum, tmp, MPFR_RNDN); 2521 mpfr_abs(tmp, tmp, MPFR_RNDN); 2522 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 2523 mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN); 2524 ++k; 2525 /* Only need to evaluate every other point on refined levels */ 2526 if (l != 1) ++k; 2527 mpfr_log10(tmp, wk, MPFR_RNDN); 2528 mpfr_abs(tmp, tmp, MPFR_RNDN); 2529 } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor * digits); /* Only need to evaluate sum until weights are < 32 digits of precision */ 2530 mpfr_sub(tmp, sum, osum, MPFR_RNDN); 2531 mpfr_abs(tmp, tmp, MPFR_RNDN); 2532 mpfr_log10(tmp, tmp, MPFR_RNDN); 2533 d1 = mpfr_get_d(tmp, MPFR_RNDN); 2534 mpfr_sub(tmp, sum, psum, MPFR_RNDN); 2535 mpfr_abs(tmp, tmp, MPFR_RNDN); 2536 mpfr_log10(tmp, tmp, MPFR_RNDN); 2537 d2 = mpfr_get_d(tmp, MPFR_RNDN); 2538 mpfr_log10(tmp, maxTerm, MPFR_RNDN); 2539 d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits; 2540 mpfr_log10(tmp, curTerm, MPFR_RNDN); 2541 d4 = mpfr_get_d(tmp, MPFR_RNDN); 2542 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4))); 2543 } while (d < digits && l < 8); 2544 *sol = mpfr_get_d(sum, MPFR_RNDN); 2545 /* Cleanup */ 2546 mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL); 2547 PetscFunctionReturn(PETSC_SUCCESS); 2548 } 2549 #else 2550 2551 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol) 2552 { 2553 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp"); 2554 } 2555 #endif 2556 2557 /*@ 2558 PetscDTTensorQuadratureCreate - create the tensor product quadrature from two lower-dimensional quadratures 2559 2560 Not Collective 2561 2562 Input Parameters: 2563 + q1 - The first quadrature 2564 - q2 - The second quadrature 2565 2566 Output Parameter: 2567 . q - A `PetscQuadrature` object 2568 2569 Level: intermediate 2570 2571 .seealso: `PetscQuadrature`, `PetscDTGaussTensorQuadrature()` 2572 @*/ 2573 PetscErrorCode PetscDTTensorQuadratureCreate(PetscQuadrature q1, PetscQuadrature q2, PetscQuadrature *q) 2574 { 2575 DMPolytopeType ct1, ct2, ct; 2576 const PetscReal *x1, *w1, *x2, *w2; 2577 PetscReal *x, *w; 2578 PetscInt dim1, Nc1, Np1, order1, qa, d1; 2579 PetscInt dim2, Nc2, Np2, order2, qb, d2; 2580 PetscInt dim, Nc, Np, order, qc, d; 2581 2582 PetscFunctionBegin; 2583 PetscValidHeaderSpecific(q1, PETSCQUADRATURE_CLASSID, 1); 2584 PetscValidHeaderSpecific(q2, PETSCQUADRATURE_CLASSID, 2); 2585 PetscValidPointer(q, 3); 2586 PetscCall(PetscQuadratureGetOrder(q1, &order1)); 2587 PetscCall(PetscQuadratureGetOrder(q2, &order2)); 2588 PetscCheck(order1 == order2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Order1 %" PetscInt_FMT " != %" PetscInt_FMT " Order2", order1, order2); 2589 PetscCall(PetscQuadratureGetData(q1, &dim1, &Nc1, &Np1, &x1, &w1)); 2590 PetscCall(PetscQuadratureGetCellType(q1, &ct1)); 2591 PetscCall(PetscQuadratureGetData(q2, &dim2, &Nc2, &Np2, &x2, &w2)); 2592 PetscCall(PetscQuadratureGetCellType(q2, &ct2)); 2593 PetscCheck(Nc1 == Nc2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "NumComp1 %" PetscInt_FMT " != %" PetscInt_FMT " NumComp2", Nc1, Nc2); 2594 2595 switch (ct1) { 2596 case DM_POLYTOPE_POINT: 2597 ct = ct2; 2598 break; 2599 case DM_POLYTOPE_SEGMENT: 2600 switch (ct2) { 2601 case DM_POLYTOPE_POINT: 2602 ct = ct1; 2603 break; 2604 case DM_POLYTOPE_SEGMENT: 2605 ct = DM_POLYTOPE_QUADRILATERAL; 2606 break; 2607 case DM_POLYTOPE_TRIANGLE: 2608 ct = DM_POLYTOPE_TRI_PRISM; 2609 break; 2610 case DM_POLYTOPE_QUADRILATERAL: 2611 ct = DM_POLYTOPE_HEXAHEDRON; 2612 break; 2613 case DM_POLYTOPE_TETRAHEDRON: 2614 ct = DM_POLYTOPE_UNKNOWN; 2615 break; 2616 case DM_POLYTOPE_HEXAHEDRON: 2617 ct = DM_POLYTOPE_UNKNOWN; 2618 break; 2619 default: 2620 ct = DM_POLYTOPE_UNKNOWN; 2621 } 2622 break; 2623 case DM_POLYTOPE_TRIANGLE: 2624 switch (ct2) { 2625 case DM_POLYTOPE_POINT: 2626 ct = ct1; 2627 break; 2628 case DM_POLYTOPE_SEGMENT: 2629 ct = DM_POLYTOPE_TRI_PRISM; 2630 break; 2631 case DM_POLYTOPE_TRIANGLE: 2632 ct = DM_POLYTOPE_UNKNOWN; 2633 break; 2634 case DM_POLYTOPE_QUADRILATERAL: 2635 ct = DM_POLYTOPE_UNKNOWN; 2636 break; 2637 case DM_POLYTOPE_TETRAHEDRON: 2638 ct = DM_POLYTOPE_UNKNOWN; 2639 break; 2640 case DM_POLYTOPE_HEXAHEDRON: 2641 ct = DM_POLYTOPE_UNKNOWN; 2642 break; 2643 default: 2644 ct = DM_POLYTOPE_UNKNOWN; 2645 } 2646 break; 2647 case DM_POLYTOPE_QUADRILATERAL: 2648 switch (ct2) { 2649 case DM_POLYTOPE_POINT: 2650 ct = ct1; 2651 break; 2652 case DM_POLYTOPE_SEGMENT: 2653 ct = DM_POLYTOPE_HEXAHEDRON; 2654 break; 2655 case DM_POLYTOPE_TRIANGLE: 2656 ct = DM_POLYTOPE_UNKNOWN; 2657 break; 2658 case DM_POLYTOPE_QUADRILATERAL: 2659 ct = DM_POLYTOPE_UNKNOWN; 2660 break; 2661 case DM_POLYTOPE_TETRAHEDRON: 2662 ct = DM_POLYTOPE_UNKNOWN; 2663 break; 2664 case DM_POLYTOPE_HEXAHEDRON: 2665 ct = DM_POLYTOPE_UNKNOWN; 2666 break; 2667 default: 2668 ct = DM_POLYTOPE_UNKNOWN; 2669 } 2670 break; 2671 case DM_POLYTOPE_TETRAHEDRON: 2672 switch (ct2) { 2673 case DM_POLYTOPE_POINT: 2674 ct = ct1; 2675 break; 2676 case DM_POLYTOPE_SEGMENT: 2677 ct = DM_POLYTOPE_UNKNOWN; 2678 break; 2679 case DM_POLYTOPE_TRIANGLE: 2680 ct = DM_POLYTOPE_UNKNOWN; 2681 break; 2682 case DM_POLYTOPE_QUADRILATERAL: 2683 ct = DM_POLYTOPE_UNKNOWN; 2684 break; 2685 case DM_POLYTOPE_TETRAHEDRON: 2686 ct = DM_POLYTOPE_UNKNOWN; 2687 break; 2688 case DM_POLYTOPE_HEXAHEDRON: 2689 ct = DM_POLYTOPE_UNKNOWN; 2690 break; 2691 default: 2692 ct = DM_POLYTOPE_UNKNOWN; 2693 } 2694 break; 2695 case DM_POLYTOPE_HEXAHEDRON: 2696 switch (ct2) { 2697 case DM_POLYTOPE_POINT: 2698 ct = ct1; 2699 break; 2700 case DM_POLYTOPE_SEGMENT: 2701 ct = DM_POLYTOPE_UNKNOWN; 2702 break; 2703 case DM_POLYTOPE_TRIANGLE: 2704 ct = DM_POLYTOPE_UNKNOWN; 2705 break; 2706 case DM_POLYTOPE_QUADRILATERAL: 2707 ct = DM_POLYTOPE_UNKNOWN; 2708 break; 2709 case DM_POLYTOPE_TETRAHEDRON: 2710 ct = DM_POLYTOPE_UNKNOWN; 2711 break; 2712 case DM_POLYTOPE_HEXAHEDRON: 2713 ct = DM_POLYTOPE_UNKNOWN; 2714 break; 2715 default: 2716 ct = DM_POLYTOPE_UNKNOWN; 2717 } 2718 break; 2719 default: 2720 ct = DM_POLYTOPE_UNKNOWN; 2721 } 2722 dim = dim1 + dim2; 2723 Nc = Nc1; 2724 Np = Np1 * Np2; 2725 order = order1; 2726 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q)); 2727 PetscCall(PetscQuadratureSetCellType(*q, ct)); 2728 PetscCall(PetscQuadratureSetOrder(*q, order)); 2729 PetscCall(PetscMalloc1(Np * dim, &x)); 2730 PetscCall(PetscMalloc1(Np, &w)); 2731 for (qa = 0, qc = 0; qa < Np1; ++qa) { 2732 for (qb = 0; qb < Np2; ++qb, ++qc) { 2733 for (d1 = 0, d = 0; d1 < dim1; ++d1, ++d) x[qc * dim + d] = x1[qa * dim1 + d1]; 2734 for (d2 = 0; d2 < dim2; ++d2, ++d) x[qc * dim + d] = x2[qb * dim2 + d2]; 2735 w[qc] = w1[qa] * w2[qb]; 2736 } 2737 } 2738 PetscCall(PetscQuadratureSetData(*q, dim, Nc, Np, x, w)); 2739 PetscFunctionReturn(PETSC_SUCCESS); 2740 } 2741 2742 /* Overwrites A. Can only handle full-rank problems with m>=n 2743 A in column-major format 2744 Ainv in row-major format 2745 tau has length m 2746 worksize must be >= max(1,n) 2747 */ 2748 static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m, PetscInt mstride, PetscInt n, PetscReal *A_in, PetscReal *Ainv_out, PetscScalar *tau, PetscInt worksize, PetscScalar *work) 2749 { 2750 PetscBLASInt M, N, K, lda, ldb, ldwork, info; 2751 PetscScalar *A, *Ainv, *R, *Q, Alpha; 2752 2753 PetscFunctionBegin; 2754 #if defined(PETSC_USE_COMPLEX) 2755 { 2756 PetscInt i, j; 2757 PetscCall(PetscMalloc2(m * n, &A, m * n, &Ainv)); 2758 for (j = 0; j < n; j++) { 2759 for (i = 0; i < m; i++) A[i + m * j] = A_in[i + mstride * j]; 2760 } 2761 mstride = m; 2762 } 2763 #else 2764 A = A_in; 2765 Ainv = Ainv_out; 2766 #endif 2767 2768 PetscCall(PetscBLASIntCast(m, &M)); 2769 PetscCall(PetscBLASIntCast(n, &N)); 2770 PetscCall(PetscBLASIntCast(mstride, &lda)); 2771 PetscCall(PetscBLASIntCast(worksize, &ldwork)); 2772 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 2773 PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info)); 2774 PetscCall(PetscFPTrapPop()); 2775 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error"); 2776 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */ 2777 2778 /* Extract an explicit representation of Q */ 2779 Q = Ainv; 2780 PetscCall(PetscArraycpy(Q, A, mstride * n)); 2781 K = N; /* full rank */ 2782 PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info)); 2783 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error"); 2784 2785 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */ 2786 Alpha = 1.0; 2787 ldb = lda; 2788 PetscCallBLAS("BLAStrsm", BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb)); 2789 /* Ainv is Q, overwritten with inverse */ 2790 2791 #if defined(PETSC_USE_COMPLEX) 2792 { 2793 PetscInt i; 2794 for (i = 0; i < m * n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]); 2795 PetscCall(PetscFree2(A, Ainv)); 2796 } 2797 #endif 2798 PetscFunctionReturn(PETSC_SUCCESS); 2799 } 2800 2801 /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */ 2802 static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval, const PetscReal *x, PetscInt ndegree, const PetscInt *degrees, PetscBool Transpose, PetscReal *B) 2803 { 2804 PetscReal *Bv; 2805 PetscInt i, j; 2806 2807 PetscFunctionBegin; 2808 PetscCall(PetscMalloc1((ninterval + 1) * ndegree, &Bv)); 2809 /* Point evaluation of L_p on all the source vertices */ 2810 PetscCall(PetscDTLegendreEval(ninterval + 1, x, ndegree, degrees, Bv, NULL, NULL)); 2811 /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */ 2812 for (i = 0; i < ninterval; i++) { 2813 for (j = 0; j < ndegree; j++) { 2814 if (Transpose) B[i + ninterval * j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j]; 2815 else B[i * ndegree + j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j]; 2816 } 2817 } 2818 PetscCall(PetscFree(Bv)); 2819 PetscFunctionReturn(PETSC_SUCCESS); 2820 } 2821 2822 /*@ 2823 PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals 2824 2825 Not Collective 2826 2827 Input Parameters: 2828 + degree - degree of reconstruction polynomial 2829 . nsource - number of source intervals 2830 . sourcex - sorted coordinates of source cell boundaries (length nsource+1) 2831 . ntarget - number of target intervals 2832 - targetx - sorted coordinates of target cell boundaries (length ntarget+1) 2833 2834 Output Parameter: 2835 . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s] 2836 2837 Level: advanced 2838 2839 .seealso: `PetscDTLegendreEval()` 2840 @*/ 2841 PetscErrorCode PetscDTReconstructPoly(PetscInt degree, PetscInt nsource, const PetscReal *sourcex, PetscInt ntarget, const PetscReal *targetx, PetscReal *R) 2842 { 2843 PetscInt i, j, k, *bdegrees, worksize; 2844 PetscReal xmin, xmax, center, hscale, *sourcey, *targety, *Bsource, *Bsinv, *Btarget; 2845 PetscScalar *tau, *work; 2846 2847 PetscFunctionBegin; 2848 PetscValidRealPointer(sourcex, 3); 2849 PetscValidRealPointer(targetx, 5); 2850 PetscValidRealPointer(R, 6); 2851 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); 2852 if (PetscDefined(USE_DEBUG)) { 2853 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]); 2854 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]); 2855 } 2856 xmin = PetscMin(sourcex[0], targetx[0]); 2857 xmax = PetscMax(sourcex[nsource], targetx[ntarget]); 2858 center = (xmin + xmax) / 2; 2859 hscale = (xmax - xmin) / 2; 2860 worksize = nsource; 2861 PetscCall(PetscMalloc4(degree + 1, &bdegrees, nsource + 1, &sourcey, nsource * (degree + 1), &Bsource, worksize, &work)); 2862 PetscCall(PetscMalloc4(nsource, &tau, nsource * (degree + 1), &Bsinv, ntarget + 1, &targety, ntarget * (degree + 1), &Btarget)); 2863 for (i = 0; i <= nsource; i++) sourcey[i] = (sourcex[i] - center) / hscale; 2864 for (i = 0; i <= degree; i++) bdegrees[i] = i + 1; 2865 PetscCall(PetscDTLegendreIntegrate(nsource, sourcey, degree + 1, bdegrees, PETSC_TRUE, Bsource)); 2866 PetscCall(PetscDTPseudoInverseQR(nsource, nsource, degree + 1, Bsource, Bsinv, tau, nsource, work)); 2867 for (i = 0; i <= ntarget; i++) targety[i] = (targetx[i] - center) / hscale; 2868 PetscCall(PetscDTLegendreIntegrate(ntarget, targety, degree + 1, bdegrees, PETSC_FALSE, Btarget)); 2869 for (i = 0; i < ntarget; i++) { 2870 PetscReal rowsum = 0; 2871 for (j = 0; j < nsource; j++) { 2872 PetscReal sum = 0; 2873 for (k = 0; k < degree + 1; k++) sum += Btarget[i * (degree + 1) + k] * Bsinv[k * nsource + j]; 2874 R[i * nsource + j] = sum; 2875 rowsum += sum; 2876 } 2877 for (j = 0; j < nsource; j++) R[i * nsource + j] /= rowsum; /* normalize each row */ 2878 } 2879 PetscCall(PetscFree4(bdegrees, sourcey, Bsource, work)); 2880 PetscCall(PetscFree4(tau, Bsinv, targety, Btarget)); 2881 PetscFunctionReturn(PETSC_SUCCESS); 2882 } 2883 2884 /*@C 2885 PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points 2886 2887 Not Collective 2888 2889 Input Parameters: 2890 + n - the number of GLL nodes 2891 . nodes - the GLL nodes 2892 . weights - the GLL weights 2893 - f - the function values at the nodes 2894 2895 Output Parameter: 2896 . in - the value of the integral 2897 2898 Level: beginner 2899 2900 .seealso: `PetscDTGaussLobattoLegendreQuadrature()` 2901 @*/ 2902 PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n, PetscReal *nodes, PetscReal *weights, const PetscReal *f, PetscReal *in) 2903 { 2904 PetscInt i; 2905 2906 PetscFunctionBegin; 2907 *in = 0.; 2908 for (i = 0; i < n; i++) *in += f[i] * f[i] * weights[i]; 2909 PetscFunctionReturn(PETSC_SUCCESS); 2910 } 2911 2912 /*@C 2913 PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element 2914 2915 Not Collective 2916 2917 Input Parameters: 2918 + n - the number of GLL nodes 2919 . nodes - the GLL nodes 2920 - weights - the GLL weights 2921 2922 Output Parameter: 2923 . A - the stiffness element 2924 2925 Level: beginner 2926 2927 Notes: 2928 Destroy this with `PetscGaussLobattoLegendreElementLaplacianDestroy()` 2929 2930 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) 2931 2932 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()` 2933 @*/ 2934 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 2935 { 2936 PetscReal **A; 2937 const PetscReal *gllnodes = nodes; 2938 const PetscInt p = n - 1; 2939 PetscReal z0, z1, z2 = -1, x, Lpj, Lpr; 2940 PetscInt i, j, nn, r; 2941 2942 PetscFunctionBegin; 2943 PetscCall(PetscMalloc1(n, &A)); 2944 PetscCall(PetscMalloc1(n * n, &A[0])); 2945 for (i = 1; i < n; i++) A[i] = A[i - 1] + n; 2946 2947 for (j = 1; j < p; j++) { 2948 x = gllnodes[j]; 2949 z0 = 1.; 2950 z1 = x; 2951 for (nn = 1; nn < p; nn++) { 2952 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 2953 z0 = z1; 2954 z1 = z2; 2955 } 2956 Lpj = z2; 2957 for (r = 1; r < p; r++) { 2958 if (r == j) { 2959 A[j][j] = 2. / (3. * (1. - gllnodes[j] * gllnodes[j]) * Lpj * Lpj); 2960 } else { 2961 x = gllnodes[r]; 2962 z0 = 1.; 2963 z1 = x; 2964 for (nn = 1; nn < p; nn++) { 2965 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 2966 z0 = z1; 2967 z1 = z2; 2968 } 2969 Lpr = z2; 2970 A[r][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * Lpr * (gllnodes[j] - gllnodes[r]) * (gllnodes[j] - gllnodes[r])); 2971 } 2972 } 2973 } 2974 for (j = 1; j < p + 1; j++) { 2975 x = gllnodes[j]; 2976 z0 = 1.; 2977 z1 = x; 2978 for (nn = 1; nn < p; nn++) { 2979 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 2980 z0 = z1; 2981 z1 = z2; 2982 } 2983 Lpj = z2; 2984 A[j][0] = 4. * PetscPowRealInt(-1., p) / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. + gllnodes[j]) * (1. + gllnodes[j])); 2985 A[0][j] = A[j][0]; 2986 } 2987 for (j = 0; j < p; j++) { 2988 x = gllnodes[j]; 2989 z0 = 1.; 2990 z1 = x; 2991 for (nn = 1; nn < p; nn++) { 2992 z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.)); 2993 z0 = z1; 2994 z1 = z2; 2995 } 2996 Lpj = z2; 2997 2998 A[p][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. - gllnodes[j]) * (1. - gllnodes[j])); 2999 A[j][p] = A[p][j]; 3000 } 3001 A[0][0] = 0.5 + (((PetscReal)p) * (((PetscReal)p) + 1.) - 2.) / 6.; 3002 A[p][p] = A[0][0]; 3003 *AA = A; 3004 PetscFunctionReturn(PETSC_SUCCESS); 3005 } 3006 3007 /*@C 3008 PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element created with `PetscGaussLobattoLegendreElementLaplacianCreate()` 3009 3010 Not Collective 3011 3012 Input Parameters: 3013 + n - the number of GLL nodes 3014 . nodes - the GLL nodes 3015 . weights - the GLL weightss 3016 - A - the stiffness element 3017 3018 Level: beginner 3019 3020 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()` 3021 @*/ 3022 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 3023 { 3024 PetscFunctionBegin; 3025 PetscCall(PetscFree((*AA)[0])); 3026 PetscCall(PetscFree(*AA)); 3027 *AA = NULL; 3028 PetscFunctionReturn(PETSC_SUCCESS); 3029 } 3030 3031 /*@C 3032 PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element 3033 3034 Not Collective 3035 3036 Input Parameter: 3037 + n - the number of GLL nodes 3038 . nodes - the GLL nodes 3039 . weights - the GLL weights 3040 3041 Output Parameters: 3042 . AA - the stiffness element 3043 - AAT - the transpose of AA (pass in `NULL` if you do not need this array) 3044 3045 Level: beginner 3046 3047 Notes: 3048 Destroy this with `PetscGaussLobattoLegendreElementGradientDestroy()` 3049 3050 You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented 3051 3052 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`, `PetscGaussLobattoLegendreElementGradientDestroy()` 3053 @*/ 3054 PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA, PetscReal ***AAT) 3055 { 3056 PetscReal **A, **AT = NULL; 3057 const PetscReal *gllnodes = nodes; 3058 const PetscInt p = n - 1; 3059 PetscReal Li, Lj, d0; 3060 PetscInt i, j; 3061 3062 PetscFunctionBegin; 3063 PetscCall(PetscMalloc1(n, &A)); 3064 PetscCall(PetscMalloc1(n * n, &A[0])); 3065 for (i = 1; i < n; i++) A[i] = A[i - 1] + n; 3066 3067 if (AAT) { 3068 PetscCall(PetscMalloc1(n, &AT)); 3069 PetscCall(PetscMalloc1(n * n, &AT[0])); 3070 for (i = 1; i < n; i++) AT[i] = AT[i - 1] + n; 3071 } 3072 3073 if (n == 1) A[0][0] = 0.; 3074 d0 = (PetscReal)p * ((PetscReal)p + 1.) / 4.; 3075 for (i = 0; i < n; i++) { 3076 for (j = 0; j < n; j++) { 3077 A[i][j] = 0.; 3078 PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li)); 3079 PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj)); 3080 if (i != j) A[i][j] = Li / (Lj * (gllnodes[i] - gllnodes[j])); 3081 if ((j == i) && (i == 0)) A[i][j] = -d0; 3082 if (j == i && i == p) A[i][j] = d0; 3083 if (AT) AT[j][i] = A[i][j]; 3084 } 3085 } 3086 if (AAT) *AAT = AT; 3087 *AA = A; 3088 PetscFunctionReturn(PETSC_SUCCESS); 3089 } 3090 3091 /*@C 3092 PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with `PetscGaussLobattoLegendreElementGradientCreate()` 3093 3094 Not Collective 3095 3096 Input Parameters: 3097 + n - the number of GLL nodes 3098 . nodes - the GLL nodes 3099 . weights - the GLL weights 3100 . AA - the stiffness element 3101 - AAT - the transpose of the element 3102 3103 Level: beginner 3104 3105 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionCreate()` 3106 @*/ 3107 PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA, PetscReal ***AAT) 3108 { 3109 PetscFunctionBegin; 3110 PetscCall(PetscFree((*AA)[0])); 3111 PetscCall(PetscFree(*AA)); 3112 *AA = NULL; 3113 if (AAT) { 3114 PetscCall(PetscFree((*AAT)[0])); 3115 PetscCall(PetscFree(*AAT)); 3116 *AAT = NULL; 3117 } 3118 PetscFunctionReturn(PETSC_SUCCESS); 3119 } 3120 3121 /*@C 3122 PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element 3123 3124 Not Collective 3125 3126 Input Parameters: 3127 + n - the number of GLL nodes 3128 . nodes - the GLL nodes 3129 - weights - the GLL weightss 3130 3131 Output Parameter: 3132 . AA - the stiffness element 3133 3134 Level: beginner 3135 3136 Notes: 3137 Destroy this with `PetscGaussLobattoLegendreElementAdvectionDestroy()` 3138 3139 This is the same as the Gradient operator multiplied by the diagonal mass matrix 3140 3141 You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented 3142 3143 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionDestroy()` 3144 @*/ 3145 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 3146 { 3147 PetscReal **D; 3148 const PetscReal *gllweights = weights; 3149 const PetscInt glln = n; 3150 PetscInt i, j; 3151 3152 PetscFunctionBegin; 3153 PetscCall(PetscGaussLobattoLegendreElementGradientCreate(n, nodes, weights, &D, NULL)); 3154 for (i = 0; i < glln; i++) { 3155 for (j = 0; j < glln; j++) D[i][j] = gllweights[i] * D[i][j]; 3156 } 3157 *AA = D; 3158 PetscFunctionReturn(PETSC_SUCCESS); 3159 } 3160 3161 /*@C 3162 PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element created with `PetscGaussLobattoLegendreElementAdvectionCreate()` 3163 3164 Not Collective 3165 3166 Input Parameters: 3167 + n - the number of GLL nodes 3168 . nodes - the GLL nodes 3169 . weights - the GLL weights 3170 - A - advection 3171 3172 Level: beginner 3173 3174 .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementAdvectionCreate()` 3175 @*/ 3176 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 3177 { 3178 PetscFunctionBegin; 3179 PetscCall(PetscFree((*AA)[0])); 3180 PetscCall(PetscFree(*AA)); 3181 *AA = NULL; 3182 PetscFunctionReturn(PETSC_SUCCESS); 3183 } 3184 3185 PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 3186 { 3187 PetscReal **A; 3188 const PetscReal *gllweights = weights; 3189 const PetscInt glln = n; 3190 PetscInt i, j; 3191 3192 PetscFunctionBegin; 3193 PetscCall(PetscMalloc1(glln, &A)); 3194 PetscCall(PetscMalloc1(glln * glln, &A[0])); 3195 for (i = 1; i < glln; i++) A[i] = A[i - 1] + glln; 3196 if (glln == 1) A[0][0] = 0.; 3197 for (i = 0; i < glln; i++) { 3198 for (j = 0; j < glln; j++) { 3199 A[i][j] = 0.; 3200 if (j == i) A[i][j] = gllweights[i]; 3201 } 3202 } 3203 *AA = A; 3204 PetscFunctionReturn(PETSC_SUCCESS); 3205 } 3206 3207 PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA) 3208 { 3209 PetscFunctionBegin; 3210 PetscCall(PetscFree((*AA)[0])); 3211 PetscCall(PetscFree(*AA)); 3212 *AA = NULL; 3213 PetscFunctionReturn(PETSC_SUCCESS); 3214 } 3215 3216 /*@ 3217 PetscDTIndexToBary - convert an index into a barycentric coordinate. 3218 3219 Input Parameters: 3220 + 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) 3221 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 3222 - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum) 3223 3224 Output Parameter: 3225 . coord - will be filled with the barycentric coordinate 3226 3227 Level: beginner 3228 3229 Note: 3230 The indices map to barycentric coordinates in lexicographic order, where the first index is the 3231 least significant and the last index is the most significant. 3232 3233 .seealso: `PetscDTBaryToIndex()` 3234 @*/ 3235 PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[]) 3236 { 3237 PetscInt c, d, s, total, subtotal, nexttotal; 3238 3239 PetscFunctionBeginHot; 3240 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 3241 PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative"); 3242 if (!len) { 3243 if (!sum && !index) PetscFunctionReturn(PETSC_SUCCESS); 3244 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 3245 } 3246 for (c = 1, total = 1; c <= len; c++) { 3247 /* total is the number of ways to have a tuple of length c with sum */ 3248 if (index < total) break; 3249 total = (total * (sum + c)) / c; 3250 } 3251 PetscCheck(c <= len, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range"); 3252 for (d = c; d < len; d++) coord[d] = 0; 3253 for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) { 3254 /* subtotal is the number of ways to have a tuple of length c with sum s */ 3255 /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */ 3256 if ((index + subtotal) >= total) { 3257 coord[--c] = sum - s; 3258 index -= (total - subtotal); 3259 sum = s; 3260 total = nexttotal; 3261 subtotal = 1; 3262 nexttotal = 1; 3263 s = 0; 3264 } else { 3265 subtotal = (subtotal * (c + s)) / (s + 1); 3266 nexttotal = (nexttotal * (c - 1 + s)) / (s + 1); 3267 s++; 3268 } 3269 } 3270 PetscFunctionReturn(PETSC_SUCCESS); 3271 } 3272 3273 /*@ 3274 PetscDTBaryToIndex - convert a barycentric coordinate to an index 3275 3276 Input Parameters: 3277 + 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) 3278 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 3279 - coord - a barycentric coordinate with the given length and sum 3280 3281 Output Parameter: 3282 . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum) 3283 3284 Level: beginner 3285 3286 Note: 3287 The indices map to barycentric coordinates in lexicographic order, where the first index is the 3288 least significant and the last index is the most significant. 3289 3290 .seealso: `PetscDTIndexToBary` 3291 @*/ 3292 PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index) 3293 { 3294 PetscInt c; 3295 PetscInt i; 3296 PetscInt total; 3297 3298 PetscFunctionBeginHot; 3299 PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 3300 if (!len) { 3301 if (!sum) { 3302 *index = 0; 3303 PetscFunctionReturn(PETSC_SUCCESS); 3304 } 3305 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 3306 } 3307 for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c; 3308 i = total - 1; 3309 c = len - 1; 3310 sum -= coord[c]; 3311 while (sum > 0) { 3312 PetscInt subtotal; 3313 PetscInt s; 3314 3315 for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s; 3316 i -= subtotal; 3317 sum -= coord[--c]; 3318 } 3319 *index = i; 3320 PetscFunctionReturn(PETSC_SUCCESS); 3321 } 3322 3323 /*@ 3324 PetscQuadratureComputePermutations - Compute permutations of quadrature points corresponding to domain orientations 3325 3326 Input Parameter: 3327 . quad - The `PetscQuadrature` 3328 3329 Output Parameters: 3330 + Np - The number of domain orientations 3331 - perm - An array of `IS` permutations, one for ech orientation, 3332 3333 Level: developer 3334 3335 .seealso: `PetscQuadratureSetCellType()`, `PetscQuadrature` 3336 @*/ 3337 PetscErrorCode PetscQuadratureComputePermutations(PetscQuadrature quad, PetscInt *Np, IS *perm[]) 3338 { 3339 DMPolytopeType ct; 3340 const PetscReal *xq, *wq; 3341 PetscInt dim, qdim, d, Na, o, Nq, q, qp; 3342 3343 PetscFunctionBegin; 3344 PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &xq, &wq)); 3345 PetscCall(PetscQuadratureGetCellType(quad, &ct)); 3346 dim = DMPolytopeTypeGetDim(ct); 3347 Na = DMPolytopeTypeGetNumArrangments(ct); 3348 PetscCall(PetscMalloc1(Na, perm)); 3349 if (Np) *Np = Na; 3350 Na /= 2; 3351 for (o = -Na; o < Na; ++o) { 3352 DM refdm; 3353 PetscInt *idx; 3354 PetscReal xi0[3] = {-1., -1., -1.}, v0[3], J[9], detJ, txq[3]; 3355 PetscBool flg; 3356 3357 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm)); 3358 PetscCall(DMPlexOrientPoint(refdm, 0, o)); 3359 PetscCall(DMPlexComputeCellGeometryFEM(refdm, 0, NULL, v0, J, NULL, &detJ)); 3360 PetscCall(PetscMalloc1(Nq, &idx)); 3361 for (q = 0; q < Nq; ++q) { 3362 CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q * dim], txq); 3363 for (qp = 0; qp < Nq; ++qp) { 3364 PetscReal diff = 0.; 3365 3366 for (d = 0; d < dim; ++d) diff += PetscAbsReal(txq[d] - xq[qp * dim + d]); 3367 if (diff < PETSC_SMALL) break; 3368 } 3369 PetscCheck(qp < Nq, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Transformed quad point %" PetscInt_FMT " does not match another quad point", q); 3370 idx[q] = qp; 3371 } 3372 PetscCall(DMDestroy(&refdm)); 3373 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nq, idx, PETSC_OWN_POINTER, &(*perm)[o + Na])); 3374 PetscCall(ISGetInfo((*perm)[o + Na], IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg)); 3375 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ordering for orientation %" PetscInt_FMT " was not a permutation", o); 3376 PetscCall(ISSetPermutation((*perm)[o + Na])); 3377 } 3378 if (!Na) (*perm)[0] = NULL; 3379 PetscFunctionReturn(PETSC_SUCCESS); 3380 } 3381 3382 /*@ 3383 PetscDTCreateDefaultQuadrature - Create default quadrature for a given cell 3384 3385 Not collective 3386 3387 Input Parameters: 3388 + ct - The integration domain 3389 - qorder - The desired quadrature order 3390 3391 Output Parameters: 3392 + q - The cell quadrature 3393 - fq - The face quadrature 3394 3395 Level: developer 3396 3397 .seealso: `PetscFECreateDefault()`, `PetscDTGaussTensorQuadrature()`, `PetscDTSimplexQuadrature()`, `PetscDTTensorQuadratureCreate()` 3398 @*/ 3399 PetscErrorCode PetscDTCreateDefaultQuadrature(DMPolytopeType ct, PetscInt qorder, PetscQuadrature *q, PetscQuadrature *fq) 3400 { 3401 const PetscInt quadPointsPerEdge = PetscMax(qorder + 1, 1); 3402 const PetscInt dim = DMPolytopeTypeGetDim(ct); 3403 3404 PetscFunctionBegin; 3405 switch (ct) { 3406 case DM_POLYTOPE_SEGMENT: 3407 case DM_POLYTOPE_POINT_PRISM_TENSOR: 3408 case DM_POLYTOPE_QUADRILATERAL: 3409 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3410 case DM_POLYTOPE_HEXAHEDRON: 3411 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3412 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q)); 3413 PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq)); 3414 break; 3415 case DM_POLYTOPE_TRIANGLE: 3416 case DM_POLYTOPE_TETRAHEDRON: 3417 PetscCall(PetscDTSimplexQuadrature(dim, 2 * qorder, PETSCDTSIMPLEXQUAD_DEFAULT, q)); 3418 PetscCall(PetscDTSimplexQuadrature(dim - 1, 2 * qorder, PETSCDTSIMPLEXQUAD_DEFAULT, fq)); 3419 break; 3420 case DM_POLYTOPE_TRI_PRISM: 3421 case DM_POLYTOPE_TRI_PRISM_TENSOR: { 3422 PetscQuadrature q1, q2; 3423 3424 // TODO: this should be able to use symmetric rules, but doing so causes tests to fail 3425 PetscCall(PetscDTSimplexQuadrature(2, 2 * qorder, PETSCDTSIMPLEXQUAD_CONIC, &q1)); 3426 PetscCall(PetscDTGaussTensorQuadrature(1, 1, quadPointsPerEdge, -1.0, 1.0, &q2)); 3427 PetscCall(PetscDTTensorQuadratureCreate(q1, q2, q)); 3428 PetscCall(PetscQuadratureDestroy(&q2)); 3429 *fq = q1; 3430 /* TODO Need separate quadratures for each face */ 3431 } break; 3432 default: 3433 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No quadrature for celltype %s", DMPolytopeTypes[PetscMin(ct, DM_POLYTOPE_UNKNOWN)]); 3434 } 3435 PetscFunctionReturn(PETSC_SUCCESS); 3436 } 3437