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