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