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 <petscviewer.h> 8 #include <petscdmplex.h> 9 #include <petscdmshell.h> 10 11 #if defined(PETSC_HAVE_MPFR) 12 #include <mpfr.h> 13 #endif 14 15 const char *const PetscDTNodeTypes[] = {"gaussjacobi", "equispaced", "tanhsinh", "PETSCDTNODES_", NULL}; 16 17 static PetscBool GolubWelschCite = PETSC_FALSE; 18 const char GolubWelschCitation[] = "@article{GolubWelsch1969,\n" 19 " author = {Golub and Welsch},\n" 20 " title = {Calculation of Quadrature Rules},\n" 21 " journal = {Math. Comp.},\n" 22 " volume = {23},\n" 23 " number = {106},\n" 24 " pages = {221--230},\n" 25 " year = {1969}\n}\n"; 26 27 /* Numerical tests in src/dm/dt/tests/ex1.c show that when computing the nodes and weights of Gauss-Jacobi 28 quadrature rules: 29 30 - in double precision, Newton's method and Golub & Welsch both work for moderate degrees (< 100), 31 - in single precision, Newton's method starts producing incorrect roots around n = 15, but 32 the weights from Golub & Welsch become a problem before then: they produces errors 33 in computing the Jacobi-polynomial Gram matrix around n = 6. 34 35 So we default to Newton's method (required fewer dependencies) */ 36 PetscBool PetscDTGaussQuadratureNewton_Internal = PETSC_TRUE; 37 38 PetscClassId PETSCQUADRATURE_CLASSID = 0; 39 40 /*@ 41 PetscQuadratureCreate - Create a PetscQuadrature object 42 43 Collective 44 45 Input Parameter: 46 . comm - The communicator for the PetscQuadrature object 47 48 Output Parameter: 49 . q - The PetscQuadrature object 50 51 Level: beginner 52 53 .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData() 54 @*/ 55 PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q) 56 { 57 PetscErrorCode ierr; 58 59 PetscFunctionBegin; 60 PetscValidPointer(q, 2); 61 ierr = DMInitializePackage();CHKERRQ(ierr); 62 ierr = PetscHeaderCreate(*q,PETSCQUADRATURE_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);CHKERRQ(ierr); 63 (*q)->dim = -1; 64 (*q)->Nc = 1; 65 (*q)->order = -1; 66 (*q)->numPoints = 0; 67 (*q)->points = NULL; 68 (*q)->weights = NULL; 69 PetscFunctionReturn(0); 70 } 71 72 /*@ 73 PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object 74 75 Collective on q 76 77 Input Parameter: 78 . q - The PetscQuadrature object 79 80 Output Parameter: 81 . r - The new PetscQuadrature object 82 83 Level: beginner 84 85 .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData() 86 @*/ 87 PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r) 88 { 89 PetscInt order, dim, Nc, Nq; 90 const PetscReal *points, *weights; 91 PetscReal *p, *w; 92 PetscErrorCode ierr; 93 94 PetscFunctionBegin; 95 PetscValidPointer(q, 1); 96 ierr = PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);CHKERRQ(ierr); 97 ierr = PetscQuadratureGetOrder(q, &order);CHKERRQ(ierr); 98 ierr = PetscQuadratureSetOrder(*r, order);CHKERRQ(ierr); 99 ierr = PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);CHKERRQ(ierr); 100 ierr = PetscMalloc1(Nq*dim, &p);CHKERRQ(ierr); 101 ierr = PetscMalloc1(Nq*Nc, &w);CHKERRQ(ierr); 102 ierr = PetscArraycpy(p, points, Nq*dim);CHKERRQ(ierr); 103 ierr = PetscArraycpy(w, weights, Nc * Nq);CHKERRQ(ierr); 104 ierr = PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);CHKERRQ(ierr); 105 PetscFunctionReturn(0); 106 } 107 108 /*@ 109 PetscQuadratureDestroy - Destroys a PetscQuadrature object 110 111 Collective on q 112 113 Input Parameter: 114 . q - The PetscQuadrature object 115 116 Level: beginner 117 118 .seealso: PetscQuadratureCreate(), PetscQuadratureGetData() 119 @*/ 120 PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q) 121 { 122 PetscErrorCode ierr; 123 124 PetscFunctionBegin; 125 if (!*q) PetscFunctionReturn(0); 126 PetscValidHeaderSpecific((*q),PETSCQUADRATURE_CLASSID,1); 127 if (--((PetscObject)(*q))->refct > 0) { 128 *q = NULL; 129 PetscFunctionReturn(0); 130 } 131 ierr = PetscFree((*q)->points);CHKERRQ(ierr); 132 ierr = PetscFree((*q)->weights);CHKERRQ(ierr); 133 ierr = PetscHeaderDestroy(q);CHKERRQ(ierr); 134 PetscFunctionReturn(0); 135 } 136 137 /*@ 138 PetscQuadratureGetOrder - Return the order of the method 139 140 Not collective 141 142 Input Parameter: 143 . q - The PetscQuadrature object 144 145 Output Parameter: 146 . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated 147 148 Level: intermediate 149 150 .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData() 151 @*/ 152 PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order) 153 { 154 PetscFunctionBegin; 155 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 156 PetscValidPointer(order, 2); 157 *order = q->order; 158 PetscFunctionReturn(0); 159 } 160 161 /*@ 162 PetscQuadratureSetOrder - Return the order of the method 163 164 Not collective 165 166 Input Parameters: 167 + q - The PetscQuadrature object 168 - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated 169 170 Level: intermediate 171 172 .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData() 173 @*/ 174 PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order) 175 { 176 PetscFunctionBegin; 177 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 178 q->order = order; 179 PetscFunctionReturn(0); 180 } 181 182 /*@ 183 PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated 184 185 Not collective 186 187 Input Parameter: 188 . q - The PetscQuadrature object 189 190 Output Parameter: 191 . Nc - The number of components 192 193 Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components. 194 195 Level: intermediate 196 197 .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData() 198 @*/ 199 PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc) 200 { 201 PetscFunctionBegin; 202 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 203 PetscValidPointer(Nc, 2); 204 *Nc = q->Nc; 205 PetscFunctionReturn(0); 206 } 207 208 /*@ 209 PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated 210 211 Not collective 212 213 Input Parameters: 214 + q - The PetscQuadrature object 215 - Nc - The number of components 216 217 Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components. 218 219 Level: intermediate 220 221 .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData() 222 @*/ 223 PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc) 224 { 225 PetscFunctionBegin; 226 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 227 q->Nc = Nc; 228 PetscFunctionReturn(0); 229 } 230 231 /*@C 232 PetscQuadratureGetData - Returns the data defining the quadrature 233 234 Not collective 235 236 Input Parameter: 237 . q - The PetscQuadrature object 238 239 Output Parameters: 240 + dim - The spatial dimension 241 . Nc - The number of components 242 . npoints - The number of quadrature points 243 . points - The coordinates of each quadrature point 244 - weights - The weight of each quadrature point 245 246 Level: intermediate 247 248 Fortran Notes: 249 From Fortran you must call PetscQuadratureRestoreData() when you are done with the data 250 251 .seealso: PetscQuadratureCreate(), PetscQuadratureSetData() 252 @*/ 253 PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[]) 254 { 255 PetscFunctionBegin; 256 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 257 if (dim) { 258 PetscValidPointer(dim, 2); 259 *dim = q->dim; 260 } 261 if (Nc) { 262 PetscValidPointer(Nc, 3); 263 *Nc = q->Nc; 264 } 265 if (npoints) { 266 PetscValidPointer(npoints, 4); 267 *npoints = q->numPoints; 268 } 269 if (points) { 270 PetscValidPointer(points, 5); 271 *points = q->points; 272 } 273 if (weights) { 274 PetscValidPointer(weights, 6); 275 *weights = q->weights; 276 } 277 PetscFunctionReturn(0); 278 } 279 280 static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[]) 281 { 282 PetscScalar *Js, *Jinvs; 283 PetscInt i, j, k; 284 PetscBLASInt bm, bn, info; 285 PetscErrorCode ierr; 286 287 PetscFunctionBegin; 288 if (!m || !n) PetscFunctionReturn(0); 289 ierr = PetscBLASIntCast(m, &bm);CHKERRQ(ierr); 290 ierr = PetscBLASIntCast(n, &bn);CHKERRQ(ierr); 291 #if defined(PETSC_USE_COMPLEX) 292 ierr = PetscMalloc2(m*n, &Js, m*n, &Jinvs);CHKERRQ(ierr); 293 for (i = 0; i < m*n; i++) Js[i] = J[i]; 294 #else 295 Js = (PetscReal *) J; 296 Jinvs = Jinv; 297 #endif 298 if (m == n) { 299 PetscBLASInt *pivots; 300 PetscScalar *W; 301 302 ierr = PetscMalloc2(m, &pivots, m, &W);CHKERRQ(ierr); 303 304 ierr = PetscArraycpy(Jinvs, Js, m * m);CHKERRQ(ierr); 305 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info)); 306 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info); 307 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &info)); 308 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info); 309 ierr = PetscFree2(pivots, W);CHKERRQ(ierr); 310 } else if (m < n) { 311 PetscScalar *JJT; 312 PetscBLASInt *pivots; 313 PetscScalar *W; 314 315 ierr = PetscMalloc1(m*m, &JJT);CHKERRQ(ierr); 316 ierr = PetscMalloc2(m, &pivots, m, &W);CHKERRQ(ierr); 317 for (i = 0; i < m; i++) { 318 for (j = 0; j < m; j++) { 319 PetscScalar val = 0.; 320 321 for (k = 0; k < n; k++) val += Js[i * n + k] * Js[j * n + k]; 322 JJT[i * m + j] = val; 323 } 324 } 325 326 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info)); 327 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info); 328 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info)); 329 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info); 330 for (i = 0; i < n; i++) { 331 for (j = 0; j < m; j++) { 332 PetscScalar val = 0.; 333 334 for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j]; 335 Jinvs[i * m + j] = val; 336 } 337 } 338 ierr = PetscFree2(pivots, W);CHKERRQ(ierr); 339 ierr = PetscFree(JJT);CHKERRQ(ierr); 340 } else { 341 PetscScalar *JTJ; 342 PetscBLASInt *pivots; 343 PetscScalar *W; 344 345 ierr = PetscMalloc1(n*n, &JTJ);CHKERRQ(ierr); 346 ierr = PetscMalloc2(n, &pivots, n, &W);CHKERRQ(ierr); 347 for (i = 0; i < n; i++) { 348 for (j = 0; j < n; j++) { 349 PetscScalar val = 0.; 350 351 for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j]; 352 JTJ[i * n + j] = val; 353 } 354 } 355 356 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bn, pivots, &info)); 357 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info); 358 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info)); 359 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info); 360 for (i = 0; i < n; i++) { 361 for (j = 0; j < m; j++) { 362 PetscScalar val = 0.; 363 364 for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k]; 365 Jinvs[i * m + j] = val; 366 } 367 } 368 ierr = PetscFree2(pivots, W);CHKERRQ(ierr); 369 ierr = PetscFree(JTJ);CHKERRQ(ierr); 370 } 371 #if defined(PETSC_USE_COMPLEX) 372 for (i = 0; i < m*n; i++) Jinv[i] = PetscRealPart(Jinvs[i]); 373 ierr = PetscFree2(Js, Jinvs);CHKERRQ(ierr); 374 #endif 375 PetscFunctionReturn(0); 376 } 377 378 /*@ 379 PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation. 380 381 Collecive on PetscQuadrature 382 383 Input Parameters: 384 + q - the quadrature functional 385 . imageDim - the dimension of the image of the transformation 386 . origin - a point in the original space 387 . originImage - the image of the origin under the transformation 388 . J - the Jacobian of the image: an [imageDim x dim] matrix in row major order 389 - 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] 390 391 Output Parameters: 392 . 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. 393 394 Note: 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. 395 396 Level: intermediate 397 398 .seealso: PetscDTAltVPullback(), PetscDTAltVPullbackMatrix() 399 @*/ 400 PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq) 401 { 402 PetscInt dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c; 403 const PetscReal *points; 404 const PetscReal *weights; 405 PetscReal *imagePoints, *imageWeights; 406 PetscReal *Jinv; 407 PetscReal *Jinvstar; 408 PetscErrorCode ierr; 409 410 PetscFunctionBegin; 411 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 412 if (imageDim < PetscAbsInt(formDegree)) SETERRQ2(PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Cannot represent a %D-form in %D dimensions", PetscAbsInt(formDegree), imageDim); 413 ierr = PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights);CHKERRQ(ierr); 414 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize);CHKERRQ(ierr); 415 if (Nc % formSize) SETERRQ2(PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Number of components %D is not a multiple of formSize %D\n", Nc, formSize); 416 Ncopies = Nc / formSize; 417 ierr = PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize);CHKERRQ(ierr); 418 imageNc = Ncopies * imageFormSize; 419 ierr = PetscMalloc1(Npoints * imageDim, &imagePoints);CHKERRQ(ierr); 420 ierr = PetscMalloc1(Npoints * imageNc, &imageWeights);CHKERRQ(ierr); 421 ierr = PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar);CHKERRQ(ierr); 422 ierr = PetscDTJacobianInverse_Internal(imageDim, dim, J, Jinv);CHKERRQ(ierr); 423 ierr = PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar);CHKERRQ(ierr); 424 for (pt = 0; pt < Npoints; pt++) { 425 const PetscReal *point = &points[pt * dim]; 426 PetscReal *imagePoint = &imagePoints[pt * imageDim]; 427 428 for (i = 0; i < imageDim; i++) { 429 PetscReal val = originImage[i]; 430 431 for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]); 432 imagePoint[i] = val; 433 } 434 for (c = 0; c < Ncopies; c++) { 435 const PetscReal *form = &weights[pt * Nc + c * formSize]; 436 PetscReal *imageForm = &imageWeights[pt * imageNc + c * imageFormSize]; 437 438 for (i = 0; i < imageFormSize; i++) { 439 PetscReal val = 0.; 440 441 for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j]; 442 imageForm[i] = val; 443 } 444 } 445 } 446 ierr = PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq);CHKERRQ(ierr); 447 ierr = PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights);CHKERRQ(ierr); 448 ierr = PetscFree2(Jinv, Jinvstar);CHKERRQ(ierr); 449 PetscFunctionReturn(0); 450 } 451 452 /*@C 453 PetscQuadratureSetData - Sets the data defining the quadrature 454 455 Not collective 456 457 Input Parameters: 458 + q - The PetscQuadrature object 459 . dim - The spatial dimension 460 . Nc - The number of components 461 . npoints - The number of quadrature points 462 . points - The coordinates of each quadrature point 463 - weights - The weight of each quadrature point 464 465 Note: This routine owns the references to points and weights, so they must be allocated using PetscMalloc() and the user should not free them. 466 467 Level: intermediate 468 469 .seealso: PetscQuadratureCreate(), PetscQuadratureGetData() 470 @*/ 471 PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[]) 472 { 473 PetscFunctionBegin; 474 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 475 if (dim >= 0) q->dim = dim; 476 if (Nc >= 0) q->Nc = Nc; 477 if (npoints >= 0) q->numPoints = npoints; 478 if (points) { 479 PetscValidPointer(points, 5); 480 q->points = points; 481 } 482 if (weights) { 483 PetscValidPointer(weights, 6); 484 q->weights = weights; 485 } 486 PetscFunctionReturn(0); 487 } 488 489 static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v) 490 { 491 PetscInt q, d, c; 492 PetscViewerFormat format; 493 PetscErrorCode ierr; 494 495 PetscFunctionBegin; 496 if (quad->Nc > 1) {ierr = PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D) with %D components\n", quad->order, quad->numPoints, quad->dim, quad->Nc);CHKERRQ(ierr);} 497 else {ierr = PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D)\n", quad->order, quad->numPoints, quad->dim);CHKERRQ(ierr);} 498 ierr = PetscViewerGetFormat(v, &format);CHKERRQ(ierr); 499 if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 500 for (q = 0; q < quad->numPoints; ++q) { 501 ierr = PetscViewerASCIIPrintf(v, "p%D (", q);CHKERRQ(ierr); 502 ierr = PetscViewerASCIIUseTabs(v, PETSC_FALSE);CHKERRQ(ierr); 503 for (d = 0; d < quad->dim; ++d) { 504 if (d) {ierr = PetscViewerASCIIPrintf(v, ", ");CHKERRQ(ierr);} 505 ierr = PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q*quad->dim+d]);CHKERRQ(ierr); 506 } 507 ierr = PetscViewerASCIIPrintf(v, ") ");CHKERRQ(ierr); 508 if (quad->Nc > 1) {ierr = PetscViewerASCIIPrintf(v, "w%D (", q);CHKERRQ(ierr);} 509 for (c = 0; c < quad->Nc; ++c) { 510 if (c) {ierr = PetscViewerASCIIPrintf(v, ", ");CHKERRQ(ierr);} 511 ierr = PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q*quad->Nc+c]);CHKERRQ(ierr); 512 } 513 if (quad->Nc > 1) {ierr = PetscViewerASCIIPrintf(v, ")");CHKERRQ(ierr);} 514 ierr = PetscViewerASCIIPrintf(v, "\n");CHKERRQ(ierr); 515 ierr = PetscViewerASCIIUseTabs(v, PETSC_TRUE);CHKERRQ(ierr); 516 } 517 PetscFunctionReturn(0); 518 } 519 520 /*@C 521 PetscQuadratureView - Views a PetscQuadrature object 522 523 Collective on quad 524 525 Input Parameters: 526 + quad - The PetscQuadrature object 527 - viewer - The PetscViewer object 528 529 Level: beginner 530 531 .seealso: PetscQuadratureCreate(), PetscQuadratureGetData() 532 @*/ 533 PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer) 534 { 535 PetscBool iascii; 536 PetscErrorCode ierr; 537 538 PetscFunctionBegin; 539 PetscValidHeader(quad, 1); 540 if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 541 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) quad), &viewer);CHKERRQ(ierr);} 542 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 543 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 544 if (iascii) {ierr = PetscQuadratureView_Ascii(quad, viewer);CHKERRQ(ierr);} 545 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 546 PetscFunctionReturn(0); 547 } 548 549 /*@C 550 PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement 551 552 Not collective 553 554 Input Parameters: 555 + q - The original PetscQuadrature 556 . numSubelements - The number of subelements the original element is divided into 557 . v0 - An array of the initial points for each subelement 558 - jac - An array of the Jacobian mappings from the reference to each subelement 559 560 Output Parameters: 561 . dim - The dimension 562 563 Note: Together v0 and jac define an affine mapping from the original reference element to each subelement 564 565 Not available from Fortran 566 567 Level: intermediate 568 569 .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension() 570 @*/ 571 PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref) 572 { 573 const PetscReal *points, *weights; 574 PetscReal *pointsRef, *weightsRef; 575 PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e; 576 PetscErrorCode ierr; 577 578 PetscFunctionBegin; 579 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 580 PetscValidPointer(v0, 3); 581 PetscValidPointer(jac, 4); 582 PetscValidPointer(qref, 5); 583 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, qref);CHKERRQ(ierr); 584 ierr = PetscQuadratureGetOrder(q, &order);CHKERRQ(ierr); 585 ierr = PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);CHKERRQ(ierr); 586 npointsRef = npoints*numSubelements; 587 ierr = PetscMalloc1(npointsRef*dim,&pointsRef);CHKERRQ(ierr); 588 ierr = PetscMalloc1(npointsRef*Nc, &weightsRef);CHKERRQ(ierr); 589 for (c = 0; c < numSubelements; ++c) { 590 for (p = 0; p < npoints; ++p) { 591 for (d = 0; d < dim; ++d) { 592 pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d]; 593 for (e = 0; e < dim; ++e) { 594 pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0); 595 } 596 } 597 /* Could also use detJ here */ 598 for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements; 599 } 600 } 601 ierr = PetscQuadratureSetOrder(*qref, order);CHKERRQ(ierr); 602 ierr = PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);CHKERRQ(ierr); 603 PetscFunctionReturn(0); 604 } 605 606 /* Compute the coefficients for the Jacobi polynomial recurrence, 607 * 608 * J^{a,b}_n(x) = (cnm1 + cnm1x * x) * J^{a,b}_{n-1}(x) - cnm2 * J^{a,b}_{n-2}(x). 609 */ 610 #define PetscDTJacobiRecurrence_Internal(n,a,b,cnm1,cnm1x,cnm2) \ 611 do { \ 612 PetscReal _a = (a); \ 613 PetscReal _b = (b); \ 614 PetscReal _n = (n); \ 615 if (n == 1) { \ 616 (cnm1) = (_a-_b) * 0.5; \ 617 (cnm1x) = (_a+_b+2.)*0.5; \ 618 (cnm2) = 0.; \ 619 } else { \ 620 PetscReal _2n = _n+_n; \ 621 PetscReal _d = (_2n*(_n+_a+_b)*(_2n+_a+_b-2)); \ 622 PetscReal _n1 = (_2n+_a+_b-1.)*(_a*_a-_b*_b); \ 623 PetscReal _n1x = (_2n+_a+_b-1.)*(_2n+_a+_b)*(_2n+_a+_b-2); \ 624 PetscReal _n2 = 2.*((_n+_a-1.)*(_n+_b-1.)*(_2n+_a+_b)); \ 625 (cnm1) = _n1 / _d; \ 626 (cnm1x) = _n1x / _d; \ 627 (cnm2) = _n2 / _d; \ 628 } \ 629 } while (0) 630 631 /*@ 632 PetscDTJacobiNorm - Compute the weighted L2 norm of a Jacobi polynomial. 633 634 $\| P^{\alpha,\beta}_n \|_{\alpha,\beta}^2 = \int_{-1}^1 (1 + x)^{\alpha} (1 - x)^{\beta} P^{\alpha,\beta}_n (x)^2 dx.$ 635 636 Input Parameters: 637 - alpha - the left exponent > -1 638 . beta - the right exponent > -1 639 + n - the polynomial degree 640 641 Output Parameter: 642 . norm - the weighted L2 norm 643 644 Level: beginner 645 646 .seealso: PetscDTJacobiEval() 647 @*/ 648 PetscErrorCode PetscDTJacobiNorm(PetscReal alpha, PetscReal beta, PetscInt n, PetscReal *norm) 649 { 650 PetscReal twoab1; 651 PetscReal gr; 652 653 PetscFunctionBegin; 654 if (alpha <= -1.) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent alpha %g <= -1. invalid\n", (double) alpha); 655 if (beta <= -1.) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent beta %g <= -1. invalid\n", (double) beta); 656 if (n < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "n %D < 0 invalid\n", n); 657 twoab1 = PetscPowReal(2., alpha + beta + 1.); 658 #if defined(PETSC_HAVE_LGAMMA) 659 if (!n) { 660 gr = PetscExpReal(PetscLGamma(alpha+1.) + PetscLGamma(beta+1.) - PetscLGamma(alpha+beta+2.)); 661 } else { 662 gr = PetscExpReal(PetscLGamma(n+alpha+1.) + PetscLGamma(n+beta+1.) - (PetscLGamma(n+1.) + PetscLGamma(n+alpha+beta+1.))) / (n+n+alpha+beta+1.); 663 } 664 #else 665 { 666 PetscInt alphai = (PetscInt) alpha; 667 PetscInt betai = (PetscInt) beta; 668 PetscInt i; 669 670 gr = n ? (1. / (n+n+alpha+beta+1.)) : 1.; 671 if ((PetscReal) alphai == alpha) { 672 if (!n) { 673 for (i = 0; i < alphai; i++) gr *= (i+1.) / (beta+i+1.); 674 gr /= (alpha+beta+1.); 675 } else { 676 for (i = 0; i < alphai; i++) gr *= (n+i+1.) / (n+beta+i+1.); 677 } 678 } else if ((PetscReal) betai == beta) { 679 if (!n) { 680 for (i = 0; i < betai; i++) gr *= (i+1.) / (alpha+i+2.); 681 gr /= (alpha+beta+1.); 682 } else { 683 for (i = 0; i < betai; i++) gr *= (n+i+1.) / (n+alpha+i+1.); 684 } 685 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable."); 686 } 687 #endif 688 *norm = PetscSqrtReal(twoab1 * gr); 689 PetscFunctionReturn(0); 690 } 691 692 static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p) 693 { 694 PetscReal ak, bk; 695 PetscReal abk1; 696 PetscInt i,l,maxdegree; 697 698 PetscFunctionBegin; 699 maxdegree = degrees[ndegree-1] - k; 700 ak = a + k; 701 bk = b + k; 702 abk1 = a + b + k + 1.; 703 if (maxdegree < 0) { 704 for (i = 0; i < npoints; i++) for (l = 0; l < ndegree; l++) p[i*ndegree+l] = 0.; 705 PetscFunctionReturn(0); 706 } 707 for (i=0; i<npoints; i++) { 708 PetscReal pm1,pm2,x; 709 PetscReal cnm1, cnm1x, cnm2; 710 PetscInt j,m; 711 712 x = points[i]; 713 pm2 = 1.; 714 PetscDTJacobiRecurrence_Internal(1,ak,bk,cnm1,cnm1x,cnm2); 715 pm1 = (cnm1 + cnm1x*x); 716 l = 0; 717 while (l < ndegree && degrees[l] - k < 0) { 718 p[l++] = 0.; 719 } 720 while (l < ndegree && degrees[l] - k == 0) { 721 p[l] = pm2; 722 for (m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5; 723 l++; 724 } 725 while (l < ndegree && degrees[l] - k == 1) { 726 p[l] = pm1; 727 for (m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5; 728 l++; 729 } 730 for (j=2; j<=maxdegree; j++) { 731 PetscReal pp; 732 733 PetscDTJacobiRecurrence_Internal(j,ak,bk,cnm1,cnm1x,cnm2); 734 pp = (cnm1 + cnm1x*x)*pm1 - cnm2*pm2; 735 pm2 = pm1; 736 pm1 = pp; 737 while (l < ndegree && degrees[l] - k == j) { 738 p[l] = pp; 739 for (m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5; 740 l++; 741 } 742 } 743 p += ndegree; 744 } 745 PetscFunctionReturn(0); 746 } 747 748 /*@ 749 PetscDTJacobiEvalJet - Evaluate the jet (function and derivatives) of the Jacobi polynomials basis up to a given degree. The Jacobi polynomials with indices $\alpha$ and $\beta$ are orthogonal with respect to the weighted inner product $\langle f, g \rangle = \int_{-1}^1 (1+x)^{\alpha} (1-x)^{\beta) f(x) g(x) dx$. 750 751 Input Parameters: 752 + alpha - the left exponent of the weight 753 . beta - the right exponetn of the weight 754 . npoints - the number of points to evaluate the polynomials at 755 . points - [npoints] array of point coordinates 756 . degree - the maximm degree polynomial space to evaluate, (degree + 1) will be evaluated total. 757 - k - the maximum derivative to evaluate in the jet, (k + 1) will be evaluated total. 758 759 Output Argments: 760 - p - an array containing the evaluations of the Jacobi polynomials's jets on the points. the size is (degree + 1) x 761 (k + 1) x npoints, which also describes the order of the dimensions of this three-dimensional array: the first 762 (slowest varying) dimension is polynomial degree; the second dimension is derivative order; the third (fastest 763 varying) dimension is the index of the evaluation point. 764 765 Level: advanced 766 767 .seealso: PetscDTJacobiEval(), PetscDTPKDEvalJet() 768 @*/ 769 PetscErrorCode PetscDTJacobiEvalJet(PetscReal alpha, PetscReal beta, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[]) 770 { 771 PetscInt i, j, l; 772 PetscInt *degrees; 773 PetscReal *psingle; 774 PetscErrorCode ierr; 775 776 PetscFunctionBegin; 777 if (degree == 0) { 778 PetscInt zero = 0; 779 780 for (i = 0; i <= k; i++) { 781 ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, 1, &zero, &p[i*npoints]);CHKERRQ(ierr); 782 } 783 PetscFunctionReturn(0); 784 } 785 ierr = PetscMalloc1(degree + 1, °rees);CHKERRQ(ierr); 786 ierr = PetscMalloc1((degree + 1) * npoints, &psingle);CHKERRQ(ierr); 787 for (i = 0; i <= degree; i++) degrees[i] = i; 788 for (i = 0; i <= k; i++) { 789 ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, degree + 1, degrees, psingle);CHKERRQ(ierr); 790 for (j = 0; j <= degree; j++) { 791 for (l = 0; l < npoints; l++) { 792 p[(j * (k + 1) + i) * npoints + l] = psingle[l * (degree + 1) + j]; 793 } 794 } 795 } 796 ierr = PetscFree(psingle);CHKERRQ(ierr); 797 ierr = PetscFree(degrees);CHKERRQ(ierr); 798 PetscFunctionReturn(0); 799 } 800 801 /*@ 802 PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$ 803 at points 804 805 Not Collective 806 807 Input Parameters: 808 + npoints - number of spatial points to evaluate at 809 . alpha - the left exponent > -1 810 . beta - the right exponent > -1 811 . points - array of locations to evaluate at 812 . ndegree - number of basis degrees to evaluate 813 - degrees - sorted array of degrees to evaluate 814 815 Output Parameters: 816 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL) 817 . D - row-oriented derivative evaluation matrix (or NULL) 818 - D2 - row-oriented second derivative evaluation matrix (or NULL) 819 820 Level: intermediate 821 822 .seealso: PetscDTGaussQuadrature() 823 @*/ 824 PetscErrorCode PetscDTJacobiEval(PetscInt npoints,PetscReal alpha, PetscReal beta, const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2) 825 { 826 PetscErrorCode ierr; 827 828 PetscFunctionBegin; 829 if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1."); 830 if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1."); 831 if (!npoints || !ndegree) PetscFunctionReturn(0); 832 if (B) {ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B);CHKERRQ(ierr);} 833 if (D) {ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D);CHKERRQ(ierr);} 834 if (D2) {ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2);CHKERRQ(ierr);} 835 PetscFunctionReturn(0); 836 } 837 838 /*@ 839 PetscDTLegendreEval - evaluate Legendre polynomials at points 840 841 Not Collective 842 843 Input Parameters: 844 + npoints - number of spatial points to evaluate at 845 . points - array of locations to evaluate at 846 . ndegree - number of basis degrees to evaluate 847 - degrees - sorted array of degrees to evaluate 848 849 Output Parameters: 850 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL) 851 . D - row-oriented derivative evaluation matrix (or NULL) 852 - D2 - row-oriented second derivative evaluation matrix (or NULL) 853 854 Level: intermediate 855 856 .seealso: PetscDTGaussQuadrature() 857 @*/ 858 PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2) 859 { 860 PetscErrorCode ierr; 861 862 PetscFunctionBegin; 863 ierr = PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2);CHKERRQ(ierr); 864 PetscFunctionReturn(0); 865 } 866 867 /*@ 868 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) 869 870 Input Parameters: 871 + len - the desired length of the degree tuple 872 - index - the index to convert: should be >= 0 873 874 Output Parameter: 875 . degtup - will be filled with a tuple of degrees 876 877 Level: beginner 878 879 Note: for two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples 880 acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the 881 last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1). 882 883 .seealso: PetscDTGradedOrderToIndex() 884 @*/ 885 PetscErrorCode PetscDTIndexToGradedOrder(PetscInt len, PetscInt index, PetscInt degtup[]) 886 { 887 PetscInt i, total; 888 PetscInt sum; 889 890 PetscFunctionBeginHot; 891 if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 892 if (index < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative"); 893 total = 1; 894 sum = 0; 895 while (index >= total) { 896 index -= total; 897 total = (total * (len + sum)) / (sum + 1); 898 sum++; 899 } 900 for (i = 0; i < len; i++) { 901 PetscInt c; 902 903 degtup[i] = sum; 904 for (c = 0, total = 1; c < sum; c++) { 905 /* going into the loop, total is the number of way to have a tuple of sum exactly c with length len - 1 - i */ 906 if (index < total) break; 907 index -= total; 908 total = (total * (len - 1 - i + c)) / (c + 1); 909 degtup[i]--; 910 } 911 sum -= degtup[i]; 912 } 913 PetscFunctionReturn(0); 914 } 915 916 /*@ 917 PetscDTGradedOrderToIndex - convert a tuple into an index in a graded order, the inverse of PetscDTIndexToGradedOrder(). 918 919 Input Parameters: 920 + len - the length of the degree tuple 921 - degtup - tuple with this length 922 923 Output Parameter: 924 . index - index in graded order: >= 0 925 926 Level: Beginner 927 928 Note: for two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples 929 acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the 930 last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1). 931 932 .seealso: PetscDTIndexToGradedOrder() 933 @*/ 934 PetscErrorCode PetscDTGradedOrderToIndex(PetscInt len, const PetscInt degtup[], PetscInt *index) 935 { 936 PetscInt i, idx, sum, total; 937 938 PetscFunctionBeginHot; 939 if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 940 for (i = 0, sum = 0; i < len; i++) sum += degtup[i]; 941 idx = 0; 942 total = 1; 943 for (i = 0; i < sum; i++) { 944 idx += total; 945 total = (total * (len + i)) / (i + 1); 946 } 947 for (i = 0; i < len - 1; i++) { 948 PetscInt c; 949 950 total = 1; 951 sum -= degtup[i]; 952 for (c = 0; c < sum; c++) { 953 idx += total; 954 total = (total * (len - 1 - i + c)) / (c + 1); 955 } 956 } 957 *index = idx; 958 PetscFunctionReturn(0); 959 } 960 961 static PetscBool PKDCite = PETSC_FALSE; 962 const char PKDCitation[] = "@article{Kirby2010,\n" 963 " title={Singularity-free evaluation of collapsed-coordinate orthogonal polynomials},\n" 964 " author={Kirby, Robert C},\n" 965 " journal={ACM Transactions on Mathematical Software (TOMS)},\n" 966 " volume={37},\n" 967 " number={1},\n" 968 " pages={1--16},\n" 969 " year={2010},\n" 970 " publisher={ACM New York, NY, USA}\n}\n"; 971 972 /*@ 973 PetscDTPKDEvalJet - Evaluate the jet (function and derivatives) of the Prioriol-Koornwinder-Dubiner (PKD) basis for 974 the space of polynomials up to a given degree. The PKD basis is L2-orthonormal on the biunit simplex (which is used 975 as the reference element for finite elements in PETSc), which makes it a stable basis to use for evaluating 976 polynomials in that domain. 977 978 Input Parameters: 979 + dim - the number of variables in the multivariate polynomials 980 . npoints - the number of points to evaluate the polynomials at 981 . points - [npoints x dim] array of point coordinates 982 . 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. 983 - k - the maximum order partial derivative to evaluate in the jet. There are (dim + k choose dim) partial derivatives 984 in the jet. Choosing k = 0 means to evaluate just the function and no derivatives 985 986 Output Argments: 987 - p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is ((dim + degree) 988 choose dim) x ((dim + k) choose dim) x npoints, which also describes the order of the dimensions of this 989 three-dimensional array: the first (slowest varying) dimension is basis function index; the second dimension is jet 990 index; the third (fastest varying) dimension is the index of the evaluation point. 991 992 Level: advanced 993 994 Note: The ordering of the basis functions, and the ordering of the derivatives in the jet, both follow the graded 995 ordering of PetscDTIndexToGradedOrder() and PetscDTGradedOrderToIndex(). For example, in 3D, the polynomial with 996 leading monomial x^3,y^1,z^2, which as degree tuple (2,0,1), which by PetscDTGradedOrderToIndex() has index 12 (it is the 13th basis function in the space); 997 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). 998 999 The implementation uses Kirby's singularity-free evaluation algorithm, https://doi.org/10.1145/1644001.1644006. 1000 1001 .seealso: PetscDTGradedOrderToIndex(), PetscDTIndexToGradedOrder(), PetscDTJacobiEvalJet() 1002 @*/ 1003 PetscErrorCode PetscDTPKDEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[]) 1004 { 1005 PetscInt degidx, kidx, d, pt; 1006 PetscInt Nk, Ndeg; 1007 PetscInt *ktup, *degtup; 1008 PetscReal *scales, initscale, scaleexp; 1009 PetscErrorCode ierr; 1010 1011 PetscFunctionBegin; 1012 ierr = PetscCitationsRegister(PKDCitation, &PKDCite);CHKERRQ(ierr); 1013 ierr = PetscDTBinomialInt(dim + k, k, &Nk);CHKERRQ(ierr); 1014 ierr = PetscDTBinomialInt(degree + dim, degree, &Ndeg);CHKERRQ(ierr); 1015 ierr = PetscMalloc2(dim, °tup, dim, &ktup);CHKERRQ(ierr); 1016 ierr = PetscMalloc1(Ndeg, &scales);CHKERRQ(ierr); 1017 initscale = 1.; 1018 if (dim > 1) { 1019 ierr = PetscDTBinomial(dim,2,&scaleexp);CHKERRQ(ierr); 1020 initscale = PetscPowReal(2.,scaleexp*0.5);CHKERRQ(ierr); 1021 } 1022 for (degidx = 0; degidx < Ndeg; degidx++) { 1023 PetscInt e, i; 1024 PetscInt m1idx = -1, m2idx = -1; 1025 PetscInt n; 1026 PetscInt degsum; 1027 PetscReal alpha; 1028 PetscReal cnm1, cnm1x, cnm2; 1029 PetscReal norm; 1030 1031 ierr = PetscDTIndexToGradedOrder(dim, degidx, degtup);CHKERRQ(ierr); 1032 for (d = dim - 1; d >= 0; d--) if (degtup[d]) break; 1033 if (d < 0) { /* constant is 1 everywhere, all derivatives are zero */ 1034 scales[degidx] = initscale; 1035 for (e = 0; e < dim; e++) { 1036 ierr = PetscDTJacobiNorm(e,0.,0,&norm);CHKERRQ(ierr); 1037 scales[degidx] /= norm; 1038 } 1039 for (i = 0; i < npoints; i++) p[degidx * Nk * npoints + i] = 1.; 1040 for (i = 0; i < (Nk - 1) * npoints; i++) p[(degidx * Nk + 1) * npoints + i] = 0.; 1041 continue; 1042 } 1043 n = degtup[d]; 1044 degtup[d]--; 1045 ierr = PetscDTGradedOrderToIndex(dim, degtup, &m1idx);CHKERRQ(ierr); 1046 if (degtup[d] > 0) { 1047 degtup[d]--; 1048 ierr = PetscDTGradedOrderToIndex(dim, degtup, &m2idx);CHKERRQ(ierr); 1049 degtup[d]++; 1050 } 1051 degtup[d]++; 1052 for (e = 0, degsum = 0; e < d; e++) degsum += degtup[e]; 1053 alpha = 2 * degsum + d; 1054 PetscDTJacobiRecurrence_Internal(n,alpha,0.,cnm1,cnm1x,cnm2); 1055 1056 scales[degidx] = initscale; 1057 for (e = 0, degsum = 0; e < dim; e++) { 1058 PetscInt f; 1059 PetscReal ealpha; 1060 PetscReal enorm; 1061 1062 ealpha = 2 * degsum + e; 1063 for (f = 0; f < degsum; f++) scales[degidx] *= 2.; 1064 ierr = PetscDTJacobiNorm(ealpha,0.,degtup[e],&enorm);CHKERRQ(ierr); 1065 scales[degidx] /= enorm; 1066 degsum += degtup[e]; 1067 } 1068 1069 for (pt = 0; pt < npoints; pt++) { 1070 /* compute the multipliers */ 1071 PetscReal thetanm1, thetanm1x, thetanm2; 1072 1073 thetanm1x = dim - (d+1) + 2.*points[pt * dim + d]; 1074 for (e = d+1; e < dim; e++) thetanm1x += points[pt * dim + e]; 1075 thetanm1x *= 0.5; 1076 thetanm1 = (2. - (dim-(d+1))); 1077 for (e = d+1; e < dim; e++) thetanm1 -= points[pt * dim + e]; 1078 thetanm1 *= 0.5; 1079 thetanm2 = thetanm1 * thetanm1; 1080 1081 for (kidx = 0; kidx < Nk; kidx++) { 1082 PetscInt f; 1083 1084 ierr = PetscDTIndexToGradedOrder(dim, kidx, ktup);CHKERRQ(ierr); 1085 /* first sum in the same derivative terms */ 1086 p[(degidx * Nk + kidx) * npoints + pt] = (cnm1 * thetanm1 + cnm1x * thetanm1x) * p[(m1idx * Nk + kidx) * npoints + pt]; 1087 if (m2idx >= 0) { 1088 p[(degidx * Nk + kidx) * npoints + pt] -= cnm2 * thetanm2 * p[(m2idx * Nk + kidx) * npoints + pt]; 1089 } 1090 1091 for (f = d; f < dim; f++) { 1092 PetscInt km1idx, mplty = ktup[f]; 1093 1094 if (!mplty) continue; 1095 ktup[f]--; 1096 ierr = PetscDTGradedOrderToIndex(dim, ktup, &km1idx);CHKERRQ(ierr); 1097 1098 /* the derivative of cnm1x * thetanm1x wrt x variable f is 0.5 * cnm1x if f > d otherwise it is cnm1x */ 1099 /* the derivative of cnm1 * thetanm1 wrt x variable f is 0 if f == d, otherwise it is -0.5 * cnm1 */ 1100 /* the derivative of -cnm2 * thetanm2 wrt x variable f is 0 if f == d, otherwise it is cnm2 * thetanm1 */ 1101 if (f > d) { 1102 PetscInt f2; 1103 1104 p[(degidx * Nk + kidx) * npoints + pt] += mplty * 0.5 * (cnm1x - cnm1) * p[(m1idx * Nk + km1idx) * npoints + pt]; 1105 if (m2idx >= 0) { 1106 p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm2 * thetanm1 * p[(m2idx * Nk + km1idx) * npoints + pt]; 1107 /* second derivatives of -cnm2 * thetanm2 wrt x variable f,f2 is like - 0.5 * cnm2 */ 1108 for (f2 = f; f2 < dim; f2++) { 1109 PetscInt km2idx, mplty2 = ktup[f2]; 1110 PetscInt factor; 1111 1112 if (!mplty2) continue; 1113 ktup[f2]--; 1114 ierr = PetscDTGradedOrderToIndex(dim, ktup, &km2idx);CHKERRQ(ierr); 1115 1116 factor = mplty * mplty2; 1117 if (f == f2) factor /= 2; 1118 p[(degidx * Nk + kidx) * npoints + pt] -= 0.5 * factor * cnm2 * p[(m2idx * Nk + km2idx) * npoints + pt]; 1119 ktup[f2]++; 1120 } 1121 } 1122 } else { 1123 p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm1x * p[(m1idx * Nk + km1idx) * npoints + pt]; 1124 } 1125 ktup[f]++; 1126 } 1127 } 1128 } 1129 } 1130 for (degidx = 0; degidx < Ndeg; degidx++) { 1131 PetscReal scale = scales[degidx]; 1132 PetscInt i; 1133 1134 for (i = 0; i < Nk * npoints; i++) p[degidx*Nk*npoints + i] *= scale; 1135 } 1136 ierr = PetscFree(scales);CHKERRQ(ierr); 1137 ierr = PetscFree2(degtup, ktup);CHKERRQ(ierr); 1138 PetscFunctionReturn(0); 1139 } 1140 1141 /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V 1142 * with lds n; diag and subdiag are overwritten */ 1143 static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[], 1144 PetscReal eigs[], PetscScalar V[]) 1145 { 1146 char jobz = 'V'; /* eigenvalues and eigenvectors */ 1147 char range = 'A'; /* all eigenvalues will be found */ 1148 PetscReal VL = 0.; /* ignored because range is 'A' */ 1149 PetscReal VU = 0.; /* ignored because range is 'A' */ 1150 PetscBLASInt IL = 0; /* ignored because range is 'A' */ 1151 PetscBLASInt IU = 0; /* ignored because range is 'A' */ 1152 PetscReal abstol = 0.; /* unused */ 1153 PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */ 1154 PetscBLASInt *isuppz; 1155 PetscBLASInt lwork, liwork; 1156 PetscReal workquery; 1157 PetscBLASInt iworkquery; 1158 PetscBLASInt *iwork; 1159 PetscBLASInt info; 1160 PetscReal *work = NULL; 1161 PetscErrorCode ierr; 1162 1163 PetscFunctionBegin; 1164 #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG) 1165 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found"); 1166 #endif 1167 ierr = PetscBLASIntCast(n, &bn);CHKERRQ(ierr); 1168 ierr = PetscBLASIntCast(n, &ldz);CHKERRQ(ierr); 1169 #if !defined(PETSC_MISSING_LAPACK_STEGR) 1170 ierr = PetscMalloc1(2 * n, &isuppz);CHKERRQ(ierr); 1171 lwork = -1; 1172 liwork = -1; 1173 PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,&workquery,&lwork,&iworkquery,&liwork,&info)); 1174 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEGR error"); 1175 lwork = (PetscBLASInt) workquery; 1176 liwork = (PetscBLASInt) iworkquery; 1177 ierr = PetscMalloc2(lwork, &work, liwork, &iwork);CHKERRQ(ierr); 1178 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1179 PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,work,&lwork,iwork,&liwork,&info)); 1180 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1181 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEGR error"); 1182 ierr = PetscFree2(work, iwork);CHKERRQ(ierr); 1183 ierr = PetscFree(isuppz);CHKERRQ(ierr); 1184 #elif !defined(PETSC_MISSING_LAPACK_STEQR) 1185 jobz = 'I'; /* Compute eigenvalues and eigenvectors of the 1186 tridiagonal matrix. Z is initialized to the identity 1187 matrix. */ 1188 ierr = PetscMalloc1(PetscMax(1,2*n-2),&work);CHKERRQ(ierr); 1189 PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&bn,diag,subdiag,V,&ldz,work,&info)); 1190 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1191 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error"); 1192 ierr = PetscFree(work);CHKERRQ(ierr); 1193 ierr = PetscArraycpy(eigs,diag,n);CHKERRQ(ierr); 1194 #endif 1195 PetscFunctionReturn(0); 1196 } 1197 1198 /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi 1199 * quadrature rules on the interval [-1, 1] */ 1200 static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw) 1201 { 1202 PetscReal twoab1; 1203 PetscInt m = n - 2; 1204 PetscReal a = alpha + 1.; 1205 PetscReal b = beta + 1.; 1206 PetscReal gra, grb; 1207 1208 PetscFunctionBegin; 1209 twoab1 = PetscPowReal(2., a + b - 1.); 1210 #if defined(PETSC_HAVE_LGAMMA) 1211 grb = PetscExpReal(2. * PetscLGamma(b+1.) + PetscLGamma(m+1.) + PetscLGamma(m+a+1.) - 1212 (PetscLGamma(m+b+1) + PetscLGamma(m+a+b+1.))); 1213 gra = PetscExpReal(2. * PetscLGamma(a+1.) + PetscLGamma(m+1.) + PetscLGamma(m+b+1.) - 1214 (PetscLGamma(m+a+1) + PetscLGamma(m+a+b+1.))); 1215 #else 1216 { 1217 PetscInt alphai = (PetscInt) alpha; 1218 PetscInt betai = (PetscInt) beta; 1219 PetscErrorCode ierr; 1220 1221 if ((PetscReal) alphai == alpha && (PetscReal) betai == beta) { 1222 PetscReal binom1, binom2; 1223 1224 ierr = PetscDTBinomial(m+b, b, &binom1);CHKERRQ(ierr); 1225 ierr = PetscDTBinomial(m+a+b, b, &binom2);CHKERRQ(ierr); 1226 grb = 1./ (binom1 * binom2); 1227 ierr = PetscDTBinomial(m+a, a, &binom1);CHKERRQ(ierr); 1228 ierr = PetscDTBinomial(m+a+b, a, &binom2);CHKERRQ(ierr); 1229 gra = 1./ (binom1 * binom2); 1230 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable."); 1231 } 1232 #endif 1233 *leftw = twoab1 * grb / b; 1234 *rightw = twoab1 * gra / a; 1235 PetscFunctionReturn(0); 1236 } 1237 1238 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x. 1239 Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */ 1240 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 1241 { 1242 PetscReal pn1, pn2; 1243 PetscReal cnm1, cnm1x, cnm2; 1244 PetscInt k; 1245 1246 PetscFunctionBegin; 1247 if (!n) {*P = 1.0; PetscFunctionReturn(0);} 1248 PetscDTJacobiRecurrence_Internal(1,a,b,cnm1,cnm1x,cnm2); 1249 pn2 = 1.; 1250 pn1 = cnm1 + cnm1x*x; 1251 if (n == 1) {*P = pn1; PetscFunctionReturn(0);} 1252 *P = 0.0; 1253 for (k = 2; k < n+1; ++k) { 1254 PetscDTJacobiRecurrence_Internal(k,a,b,cnm1,cnm1x,cnm2); 1255 1256 *P = (cnm1 + cnm1x*x)*pn1 - cnm2*pn2; 1257 pn2 = pn1; 1258 pn1 = *P; 1259 } 1260 PetscFunctionReturn(0); 1261 } 1262 1263 /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */ 1264 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P) 1265 { 1266 PetscReal nP; 1267 PetscInt i; 1268 PetscErrorCode ierr; 1269 1270 PetscFunctionBegin; 1271 *P = 0.0; 1272 if (k > n) PetscFunctionReturn(0); 1273 ierr = PetscDTComputeJacobi(a+k, b+k, n-k, x, &nP);CHKERRQ(ierr); 1274 for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5; 1275 *P = nP; 1276 PetscFunctionReturn(0); 1277 } 1278 1279 static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[]) 1280 { 1281 PetscInt maxIter = 100; 1282 PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON)); 1283 PetscReal a1, a6, gf; 1284 PetscInt k; 1285 PetscErrorCode ierr; 1286 1287 PetscFunctionBegin; 1288 1289 a1 = PetscPowReal(2.0, a+b+1); 1290 #if defined(PETSC_HAVE_LGAMMA) 1291 { 1292 PetscReal a2, a3, a4, a5; 1293 a2 = PetscLGamma(a + npoints + 1); 1294 a3 = PetscLGamma(b + npoints + 1); 1295 a4 = PetscLGamma(a + b + npoints + 1); 1296 a5 = PetscLGamma(npoints + 1); 1297 gf = PetscExpReal(a2 + a3 - (a4 + a5)); 1298 } 1299 #else 1300 { 1301 PetscInt ia, ib; 1302 1303 ia = (PetscInt) a; 1304 ib = (PetscInt) b; 1305 gf = 1.; 1306 if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */ 1307 for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k); 1308 } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */ 1309 for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k); 1310 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable."); 1311 } 1312 #endif 1313 1314 a6 = a1 * gf; 1315 /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses. 1316 Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */ 1317 for (k = 0; k < npoints; ++k) { 1318 PetscReal r = PetscCosReal(PETSC_PI * (1. - (4.*k + 3. + 2.*b) / (4.*npoints + 2.*(a + b + 1.)))), dP; 1319 PetscInt j; 1320 1321 if (k > 0) r = 0.5 * (r + x[k-1]); 1322 for (j = 0; j < maxIter; ++j) { 1323 PetscReal s = 0.0, delta, f, fp; 1324 PetscInt i; 1325 1326 for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]); 1327 ierr = PetscDTComputeJacobi(a, b, npoints, r, &f);CHKERRQ(ierr); 1328 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp);CHKERRQ(ierr); 1329 delta = f / (fp - f * s); 1330 r = r - delta; 1331 if (PetscAbsReal(delta) < eps) break; 1332 } 1333 x[k] = r; 1334 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP);CHKERRQ(ierr); 1335 w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP); 1336 } 1337 PetscFunctionReturn(0); 1338 } 1339 1340 /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi 1341 * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */ 1342 static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s) 1343 { 1344 PetscInt i; 1345 1346 PetscFunctionBegin; 1347 for (i = 0; i < nPoints; i++) { 1348 PetscReal A, B, C; 1349 1350 PetscDTJacobiRecurrence_Internal(i+1,a,b,A,B,C); 1351 d[i] = -A / B; 1352 if (i) s[i-1] *= C / B; 1353 if (i < nPoints - 1) s[i] = 1. / B; 1354 } 1355 PetscFunctionReturn(0); 1356 } 1357 1358 static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w) 1359 { 1360 PetscReal mu0; 1361 PetscReal ga, gb, gab; 1362 PetscInt i; 1363 PetscErrorCode ierr; 1364 1365 PetscFunctionBegin; 1366 ierr = PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite);CHKERRQ(ierr); 1367 1368 #if defined(PETSC_HAVE_TGAMMA) 1369 ga = PetscTGamma(a + 1); 1370 gb = PetscTGamma(b + 1); 1371 gab = PetscTGamma(a + b + 2); 1372 #else 1373 { 1374 PetscInt ia, ib; 1375 1376 ia = (PetscInt) a; 1377 ib = (PetscInt) b; 1378 if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */ 1379 ierr = PetscDTFactorial(ia, &ga);CHKERRQ(ierr); 1380 ierr = PetscDTFactorial(ib, &gb);CHKERRQ(ierr); 1381 ierr = PetscDTFactorial(ia + ib + 1, &gb);CHKERRQ(ierr); 1382 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable."); 1383 } 1384 #endif 1385 mu0 = PetscPowReal(2.,a + b + 1.) * ga * gb / gab; 1386 1387 #if defined(PETSCDTGAUSSIANQUADRATURE_EIG) 1388 { 1389 PetscReal *diag, *subdiag; 1390 PetscScalar *V; 1391 1392 ierr = PetscMalloc2(npoints, &diag, npoints, &subdiag);CHKERRQ(ierr); 1393 ierr = PetscMalloc1(npoints*npoints, &V);CHKERRQ(ierr); 1394 ierr = PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag);CHKERRQ(ierr); 1395 for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]); 1396 ierr = PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V);CHKERRQ(ierr); 1397 for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0; 1398 ierr = PetscFree(V);CHKERRQ(ierr); 1399 ierr = PetscFree2(diag, subdiag);CHKERRQ(ierr); 1400 } 1401 #else 1402 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found"); 1403 #endif 1404 { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the 1405 eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that 1406 the eigenvalues are sorted */ 1407 PetscBool sorted; 1408 1409 ierr = PetscSortedReal(npoints, x, &sorted);CHKERRQ(ierr); 1410 if (!sorted) { 1411 PetscInt *order, i; 1412 PetscReal *tmp; 1413 1414 ierr = PetscMalloc2(npoints, &order, npoints, &tmp);CHKERRQ(ierr); 1415 for (i = 0; i < npoints; i++) order[i] = i; 1416 ierr = PetscSortRealWithPermutation(npoints, x, order);CHKERRQ(ierr); 1417 ierr = PetscArraycpy(tmp, x, npoints);CHKERRQ(ierr); 1418 for (i = 0; i < npoints; i++) x[i] = tmp[order[i]]; 1419 ierr = PetscArraycpy(tmp, w, npoints);CHKERRQ(ierr); 1420 for (i = 0; i < npoints; i++) w[i] = tmp[order[i]]; 1421 ierr = PetscFree2(order, tmp);CHKERRQ(ierr); 1422 } 1423 } 1424 PetscFunctionReturn(0); 1425 } 1426 1427 static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1428 { 1429 PetscErrorCode ierr; 1430 1431 PetscFunctionBegin; 1432 if (npoints < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be positive"); 1433 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1434 if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1."); 1435 if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1."); 1436 1437 if (newton) { 1438 ierr = PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w);CHKERRQ(ierr); 1439 } else { 1440 ierr = PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w);CHKERRQ(ierr); 1441 } 1442 if (alpha == beta) { /* symmetrize */ 1443 PetscInt i; 1444 for (i = 0; i < (npoints + 1) / 2; i++) { 1445 PetscInt j = npoints - 1 - i; 1446 PetscReal xi = x[i]; 1447 PetscReal xj = x[j]; 1448 PetscReal wi = w[i]; 1449 PetscReal wj = w[j]; 1450 1451 x[i] = (xi - xj) / 2.; 1452 x[j] = (xj - xi) / 2.; 1453 w[i] = w[j] = (wi + wj) / 2.; 1454 } 1455 } 1456 PetscFunctionReturn(0); 1457 } 1458 1459 /*@ 1460 PetscDTGaussJacobiQuadrature - quadrature for the interval [a, b] with the weight function 1461 $(x-a)^\alpha (x-b)^\beta$. 1462 1463 Not collective 1464 1465 Input Parameters: 1466 + npoints - the number of points in the quadrature rule 1467 . a - the left endpoint of the interval 1468 . b - the right endpoint of the interval 1469 . alpha - the left exponent 1470 - beta - the right exponent 1471 1472 Output Parameters: 1473 + x - array of length npoints, the locations of the quadrature points 1474 - w - array of length npoints, the weights of the quadrature points 1475 1476 Level: intermediate 1477 1478 Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 1. 1479 @*/ 1480 PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1481 { 1482 PetscInt i; 1483 PetscErrorCode ierr; 1484 1485 PetscFunctionBegin; 1486 ierr = PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);CHKERRQ(ierr); 1487 if (a != -1. || b != 1.) { /* shift */ 1488 for (i = 0; i < npoints; i++) { 1489 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1490 w[i] *= (b - a) / 2.; 1491 } 1492 } 1493 PetscFunctionReturn(0); 1494 } 1495 1496 static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1497 { 1498 PetscInt i; 1499 PetscErrorCode ierr; 1500 1501 PetscFunctionBegin; 1502 if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be positive"); 1503 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1504 if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1."); 1505 if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1."); 1506 1507 x[0] = -1.; 1508 x[npoints-1] = 1.; 1509 if (npoints > 2) { 1510 ierr = PetscDTGaussJacobiQuadrature_Internal(npoints-2, alpha+1., beta+1., &x[1], &w[1], newton);CHKERRQ(ierr); 1511 } 1512 for (i = 1; i < npoints - 1; i++) { 1513 w[i] /= (1. - x[i]*x[i]); 1514 } 1515 ierr = PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints-1]);CHKERRQ(ierr); 1516 PetscFunctionReturn(0); 1517 } 1518 1519 /*@ 1520 PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval [a, b] with the weight function 1521 $(x-a)^\alpha (x-b)^\beta$, with endpoints a and b included as quadrature points. 1522 1523 Not collective 1524 1525 Input Parameters: 1526 + npoints - the number of points in the quadrature rule 1527 . a - the left endpoint of the interval 1528 . b - the right endpoint of the interval 1529 . alpha - the left exponent 1530 - beta - the right exponent 1531 1532 Output Parameters: 1533 + x - array of length npoints, the locations of the quadrature points 1534 - w - array of length npoints, the weights of the quadrature points 1535 1536 Level: intermediate 1537 1538 Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 3. 1539 @*/ 1540 PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1541 { 1542 PetscInt i; 1543 PetscErrorCode ierr; 1544 1545 PetscFunctionBegin; 1546 ierr = PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);CHKERRQ(ierr); 1547 if (a != -1. || b != 1.) { /* shift */ 1548 for (i = 0; i < npoints; i++) { 1549 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1550 w[i] *= (b - a) / 2.; 1551 } 1552 } 1553 PetscFunctionReturn(0); 1554 } 1555 1556 /*@ 1557 PetscDTGaussQuadrature - create Gauss-Legendre quadrature 1558 1559 Not Collective 1560 1561 Input Parameters: 1562 + npoints - number of points 1563 . a - left end of interval (often-1) 1564 - b - right end of interval (often +1) 1565 1566 Output Parameters: 1567 + x - quadrature points 1568 - w - quadrature weights 1569 1570 Level: intermediate 1571 1572 References: 1573 . 1. - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969. 1574 1575 .seealso: PetscDTLegendreEval() 1576 @*/ 1577 PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w) 1578 { 1579 PetscInt i; 1580 PetscErrorCode ierr; 1581 1582 PetscFunctionBegin; 1583 ierr = PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal);CHKERRQ(ierr); 1584 if (a != -1. || b != 1.) { /* shift */ 1585 for (i = 0; i < npoints; i++) { 1586 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1587 w[i] *= (b - a) / 2.; 1588 } 1589 } 1590 PetscFunctionReturn(0); 1591 } 1592 1593 /*@C 1594 PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre 1595 nodes of a given size on the domain [-1,1] 1596 1597 Not Collective 1598 1599 Input Parameters: 1600 + n - number of grid nodes 1601 - type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON 1602 1603 Output Parameters: 1604 + x - quadrature points 1605 - w - quadrature weights 1606 1607 Notes: 1608 For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not 1609 close enough to the desired solution 1610 1611 These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes 1612 1613 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 1614 1615 Level: intermediate 1616 1617 .seealso: PetscDTGaussQuadrature() 1618 1619 @*/ 1620 PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal *x,PetscReal *w) 1621 { 1622 PetscBool newton; 1623 PetscErrorCode ierr; 1624 1625 PetscFunctionBegin; 1626 if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide at least 2 grid points per element"); 1627 newton = (PetscBool) (type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON); 1628 ierr = PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton);CHKERRQ(ierr); 1629 PetscFunctionReturn(0); 1630 } 1631 1632 /*@ 1633 PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature 1634 1635 Not Collective 1636 1637 Input Parameters: 1638 + dim - The spatial dimension 1639 . Nc - The number of components 1640 . npoints - number of points in one dimension 1641 . a - left end of interval (often-1) 1642 - b - right end of interval (often +1) 1643 1644 Output Parameter: 1645 . q - A PetscQuadrature object 1646 1647 Level: intermediate 1648 1649 .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval() 1650 @*/ 1651 PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 1652 { 1653 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c; 1654 PetscReal *x, *w, *xw, *ww; 1655 PetscErrorCode ierr; 1656 1657 PetscFunctionBegin; 1658 ierr = PetscMalloc1(totpoints*dim,&x);CHKERRQ(ierr); 1659 ierr = PetscMalloc1(totpoints*Nc,&w);CHKERRQ(ierr); 1660 /* Set up the Golub-Welsch system */ 1661 switch (dim) { 1662 case 0: 1663 ierr = PetscFree(x);CHKERRQ(ierr); 1664 ierr = PetscFree(w);CHKERRQ(ierr); 1665 ierr = PetscMalloc1(1, &x);CHKERRQ(ierr); 1666 ierr = PetscMalloc1(Nc, &w);CHKERRQ(ierr); 1667 x[0] = 0.0; 1668 for (c = 0; c < Nc; ++c) w[c] = 1.0; 1669 break; 1670 case 1: 1671 ierr = PetscMalloc1(npoints,&ww);CHKERRQ(ierr); 1672 ierr = PetscDTGaussQuadrature(npoints, a, b, x, ww);CHKERRQ(ierr); 1673 for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i]; 1674 ierr = PetscFree(ww);CHKERRQ(ierr); 1675 break; 1676 case 2: 1677 ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr); 1678 ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr); 1679 for (i = 0; i < npoints; ++i) { 1680 for (j = 0; j < npoints; ++j) { 1681 x[(i*npoints+j)*dim+0] = xw[i]; 1682 x[(i*npoints+j)*dim+1] = xw[j]; 1683 for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j]; 1684 } 1685 } 1686 ierr = PetscFree2(xw,ww);CHKERRQ(ierr); 1687 break; 1688 case 3: 1689 ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr); 1690 ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr); 1691 for (i = 0; i < npoints; ++i) { 1692 for (j = 0; j < npoints; ++j) { 1693 for (k = 0; k < npoints; ++k) { 1694 x[((i*npoints+j)*npoints+k)*dim+0] = xw[i]; 1695 x[((i*npoints+j)*npoints+k)*dim+1] = xw[j]; 1696 x[((i*npoints+j)*npoints+k)*dim+2] = xw[k]; 1697 for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k]; 1698 } 1699 } 1700 } 1701 ierr = PetscFree2(xw,ww);CHKERRQ(ierr); 1702 break; 1703 default: 1704 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim); 1705 } 1706 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 1707 ierr = PetscQuadratureSetOrder(*q, 2*npoints-1);CHKERRQ(ierr); 1708 ierr = PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);CHKERRQ(ierr); 1709 ierr = PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");CHKERRQ(ierr); 1710 PetscFunctionReturn(0); 1711 } 1712 1713 /*@ 1714 PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex 1715 1716 Not Collective 1717 1718 Input Parameters: 1719 + dim - The simplex dimension 1720 . Nc - The number of components 1721 . npoints - The number of points in one dimension 1722 . a - left end of interval (often-1) 1723 - b - right end of interval (often +1) 1724 1725 Output Parameter: 1726 . q - A PetscQuadrature object 1727 1728 Level: intermediate 1729 1730 References: 1731 . 1. - Karniadakis and Sherwin. FIAT 1732 1733 Note: For dim == 1, this is Gauss-Legendre quadrature 1734 1735 .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature() 1736 @*/ 1737 PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 1738 { 1739 PetscInt totprev, totrem; 1740 PetscInt totpoints; 1741 PetscReal *p1, *w1; 1742 PetscReal *x, *w; 1743 PetscInt i, j, k, l, m, pt, c; 1744 PetscErrorCode ierr; 1745 1746 PetscFunctionBegin; 1747 if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now"); 1748 totpoints = 1; 1749 for (i = 0, totpoints = 1; i < dim; i++) totpoints *= npoints; 1750 ierr = PetscMalloc1(totpoints*dim, &x);CHKERRQ(ierr); 1751 ierr = PetscMalloc1(totpoints*Nc, &w);CHKERRQ(ierr); 1752 ierr = PetscMalloc2(npoints, &p1, npoints, &w1);CHKERRQ(ierr); 1753 for (i = 0; i < totpoints*Nc; i++) w[i] = 1.; 1754 for (i = 0, totprev = 1, totrem = totpoints / npoints; i < dim; i++) { 1755 PetscReal mul; 1756 1757 mul = PetscPowReal(2.,-i); 1758 ierr = PetscDTGaussJacobiQuadrature(npoints, -1., 1., i, 0.0, p1, w1);CHKERRQ(ierr); 1759 for (pt = 0, l = 0; l < totprev; l++) { 1760 for (j = 0; j < npoints; j++) { 1761 for (m = 0; m < totrem; m++, pt++) { 1762 for (k = 0; k < i; k++) x[pt*dim+k] = (x[pt*dim+k]+1.)*(1.-p1[j])*0.5 - 1.; 1763 x[pt * dim + i] = p1[j]; 1764 for (c = 0; c < Nc; c++) w[pt*Nc + c] *= mul * w1[j]; 1765 } 1766 } 1767 } 1768 totprev *= npoints; 1769 totrem /= npoints; 1770 } 1771 ierr = PetscFree2(p1, w1);CHKERRQ(ierr); 1772 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 1773 ierr = PetscQuadratureSetOrder(*q, 2*npoints-1);CHKERRQ(ierr); 1774 ierr = PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);CHKERRQ(ierr); 1775 ierr = PetscObjectChangeTypeName((PetscObject)*q,"StroudConical");CHKERRQ(ierr); 1776 PetscFunctionReturn(0); 1777 } 1778 1779 /*@ 1780 PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell 1781 1782 Not Collective 1783 1784 Input Parameters: 1785 + dim - The cell dimension 1786 . level - The number of points in one dimension, 2^l 1787 . a - left end of interval (often-1) 1788 - b - right end of interval (often +1) 1789 1790 Output Parameter: 1791 . q - A PetscQuadrature object 1792 1793 Level: intermediate 1794 1795 .seealso: PetscDTGaussTensorQuadrature() 1796 @*/ 1797 PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q) 1798 { 1799 const PetscInt p = 16; /* Digits of precision in the evaluation */ 1800 const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */ 1801 const PetscReal beta = (b+a)/2.; /* Center of the integration interval */ 1802 const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */ 1803 PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */ 1804 PetscReal wk = 0.5*PETSC_PI; /* Quadrature weight at x_k */ 1805 PetscReal *x, *w; 1806 PetscInt K, k, npoints; 1807 PetscErrorCode ierr; 1808 1809 PetscFunctionBegin; 1810 if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim); 1811 if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits"); 1812 /* Find K such that the weights are < 32 digits of precision */ 1813 for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) { 1814 wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h))); 1815 } 1816 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 1817 ierr = PetscQuadratureSetOrder(*q, 2*K+1);CHKERRQ(ierr); 1818 npoints = 2*K-1; 1819 ierr = PetscMalloc1(npoints*dim, &x);CHKERRQ(ierr); 1820 ierr = PetscMalloc1(npoints, &w);CHKERRQ(ierr); 1821 /* Center term */ 1822 x[0] = beta; 1823 w[0] = 0.5*alpha*PETSC_PI; 1824 for (k = 1; k < K; ++k) { 1825 wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1826 xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h)); 1827 x[2*k-1] = -alpha*xk+beta; 1828 w[2*k-1] = wk; 1829 x[2*k+0] = alpha*xk+beta; 1830 w[2*k+0] = wk; 1831 } 1832 ierr = PetscQuadratureSetData(*q, dim, 1, npoints, x, w);CHKERRQ(ierr); 1833 PetscFunctionReturn(0); 1834 } 1835 1836 PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 1837 { 1838 const PetscInt p = 16; /* Digits of precision in the evaluation */ 1839 const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */ 1840 const PetscReal beta = (b+a)/2.; /* Center of the integration interval */ 1841 PetscReal h = 1.0; /* Step size, length between x_k */ 1842 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 1843 PetscReal osum = 0.0; /* Integral on last level */ 1844 PetscReal psum = 0.0; /* Integral on the level before the last level */ 1845 PetscReal sum; /* Integral on current level */ 1846 PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 1847 PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 1848 PetscReal wk; /* Quadrature weight at x_k */ 1849 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 1850 PetscInt d; /* Digits of precision in the integral */ 1851 1852 PetscFunctionBegin; 1853 if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 1854 /* Center term */ 1855 func(beta, &lval); 1856 sum = 0.5*alpha*PETSC_PI*lval; 1857 /* */ 1858 do { 1859 PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4; 1860 PetscInt k = 1; 1861 1862 ++l; 1863 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */ 1864 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 1865 psum = osum; 1866 osum = sum; 1867 h *= 0.5; 1868 sum *= 0.5; 1869 do { 1870 wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1871 yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1872 lx = -alpha*(1.0 - yk)+beta; 1873 rx = alpha*(1.0 - yk)+beta; 1874 func(lx, &lval); 1875 func(rx, &rval); 1876 lterm = alpha*wk*lval; 1877 maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm); 1878 sum += lterm; 1879 rterm = alpha*wk*rval; 1880 maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm); 1881 sum += rterm; 1882 ++k; 1883 /* Only need to evaluate every other point on refined levels */ 1884 if (l != 1) ++k; 1885 } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */ 1886 1887 d1 = PetscLog10Real(PetscAbsReal(sum - osum)); 1888 d2 = PetscLog10Real(PetscAbsReal(sum - psum)); 1889 d3 = PetscLog10Real(maxTerm) - p; 1890 if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0; 1891 else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm))); 1892 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4))); 1893 } while (d < digits && l < 12); 1894 *sol = sum; 1895 1896 PetscFunctionReturn(0); 1897 } 1898 1899 #if defined(PETSC_HAVE_MPFR) 1900 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 1901 { 1902 const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */ 1903 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 1904 mpfr_t alpha; /* Half-width of the integration interval */ 1905 mpfr_t beta; /* Center of the integration interval */ 1906 mpfr_t h; /* Step size, length between x_k */ 1907 mpfr_t osum; /* Integral on last level */ 1908 mpfr_t psum; /* Integral on the level before the last level */ 1909 mpfr_t sum; /* Integral on current level */ 1910 mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 1911 mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 1912 mpfr_t wk; /* Quadrature weight at x_k */ 1913 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 1914 PetscInt d; /* Digits of precision in the integral */ 1915 mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp; 1916 1917 PetscFunctionBegin; 1918 if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 1919 /* Create high precision storage */ 1920 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); 1921 /* Initialization */ 1922 mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN); 1923 mpfr_set_d(beta, 0.5*(b+a), MPFR_RNDN); 1924 mpfr_set_d(osum, 0.0, MPFR_RNDN); 1925 mpfr_set_d(psum, 0.0, MPFR_RNDN); 1926 mpfr_set_d(h, 1.0, MPFR_RNDN); 1927 mpfr_const_pi(pi2, MPFR_RNDN); 1928 mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN); 1929 /* Center term */ 1930 func(0.5*(b+a), &lval); 1931 mpfr_set(sum, pi2, MPFR_RNDN); 1932 mpfr_mul(sum, sum, alpha, MPFR_RNDN); 1933 mpfr_mul_d(sum, sum, lval, MPFR_RNDN); 1934 /* */ 1935 do { 1936 PetscReal d1, d2, d3, d4; 1937 PetscInt k = 1; 1938 1939 ++l; 1940 mpfr_set_d(maxTerm, 0.0, MPFR_RNDN); 1941 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */ 1942 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 1943 mpfr_set(psum, osum, MPFR_RNDN); 1944 mpfr_set(osum, sum, MPFR_RNDN); 1945 mpfr_mul_d(h, h, 0.5, MPFR_RNDN); 1946 mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN); 1947 do { 1948 mpfr_set_si(kh, k, MPFR_RNDN); 1949 mpfr_mul(kh, kh, h, MPFR_RNDN); 1950 /* Weight */ 1951 mpfr_set(wk, h, MPFR_RNDN); 1952 mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN); 1953 mpfr_mul(msinh, msinh, pi2, MPFR_RNDN); 1954 mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN); 1955 mpfr_cosh(tmp, msinh, MPFR_RNDN); 1956 mpfr_sqr(tmp, tmp, MPFR_RNDN); 1957 mpfr_mul(wk, wk, mcosh, MPFR_RNDN); 1958 mpfr_div(wk, wk, tmp, MPFR_RNDN); 1959 /* Abscissa */ 1960 mpfr_set_d(yk, 1.0, MPFR_RNDZ); 1961 mpfr_cosh(tmp, msinh, MPFR_RNDN); 1962 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 1963 mpfr_exp(tmp, msinh, MPFR_RNDN); 1964 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 1965 /* Quadrature points */ 1966 mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ); 1967 mpfr_mul(lx, lx, alpha, MPFR_RNDU); 1968 mpfr_add(lx, lx, beta, MPFR_RNDU); 1969 mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ); 1970 mpfr_mul(rx, rx, alpha, MPFR_RNDD); 1971 mpfr_add(rx, rx, beta, MPFR_RNDD); 1972 /* Evaluation */ 1973 func(mpfr_get_d(lx, MPFR_RNDU), &lval); 1974 func(mpfr_get_d(rx, MPFR_RNDD), &rval); 1975 /* Update */ 1976 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 1977 mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN); 1978 mpfr_add(sum, sum, tmp, MPFR_RNDN); 1979 mpfr_abs(tmp, tmp, MPFR_RNDN); 1980 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 1981 mpfr_set(curTerm, tmp, MPFR_RNDN); 1982 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 1983 mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN); 1984 mpfr_add(sum, sum, tmp, MPFR_RNDN); 1985 mpfr_abs(tmp, tmp, MPFR_RNDN); 1986 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 1987 mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN); 1988 ++k; 1989 /* Only need to evaluate every other point on refined levels */ 1990 if (l != 1) ++k; 1991 mpfr_log10(tmp, wk, MPFR_RNDN); 1992 mpfr_abs(tmp, tmp, MPFR_RNDN); 1993 } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */ 1994 mpfr_sub(tmp, sum, osum, MPFR_RNDN); 1995 mpfr_abs(tmp, tmp, MPFR_RNDN); 1996 mpfr_log10(tmp, tmp, MPFR_RNDN); 1997 d1 = mpfr_get_d(tmp, MPFR_RNDN); 1998 mpfr_sub(tmp, sum, psum, MPFR_RNDN); 1999 mpfr_abs(tmp, tmp, MPFR_RNDN); 2000 mpfr_log10(tmp, tmp, MPFR_RNDN); 2001 d2 = mpfr_get_d(tmp, MPFR_RNDN); 2002 mpfr_log10(tmp, maxTerm, MPFR_RNDN); 2003 d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits; 2004 mpfr_log10(tmp, curTerm, MPFR_RNDN); 2005 d4 = mpfr_get_d(tmp, MPFR_RNDN); 2006 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4))); 2007 } while (d < digits && l < 8); 2008 *sol = mpfr_get_d(sum, MPFR_RNDN); 2009 /* Cleanup */ 2010 mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL); 2011 PetscFunctionReturn(0); 2012 } 2013 #else 2014 2015 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 2016 { 2017 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp"); 2018 } 2019 #endif 2020 2021 /* Overwrites A. Can only handle full-rank problems with m>=n 2022 * A in column-major format 2023 * Ainv in row-major format 2024 * tau has length m 2025 * worksize must be >= max(1,n) 2026 */ 2027 static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work) 2028 { 2029 PetscErrorCode ierr; 2030 PetscBLASInt M,N,K,lda,ldb,ldwork,info; 2031 PetscScalar *A,*Ainv,*R,*Q,Alpha; 2032 2033 PetscFunctionBegin; 2034 #if defined(PETSC_USE_COMPLEX) 2035 { 2036 PetscInt i,j; 2037 ierr = PetscMalloc2(m*n,&A,m*n,&Ainv);CHKERRQ(ierr); 2038 for (j=0; j<n; j++) { 2039 for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j]; 2040 } 2041 mstride = m; 2042 } 2043 #else 2044 A = A_in; 2045 Ainv = Ainv_out; 2046 #endif 2047 2048 ierr = PetscBLASIntCast(m,&M);CHKERRQ(ierr); 2049 ierr = PetscBLASIntCast(n,&N);CHKERRQ(ierr); 2050 ierr = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr); 2051 ierr = PetscBLASIntCast(worksize,&ldwork);CHKERRQ(ierr); 2052 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 2053 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info)); 2054 ierr = PetscFPTrapPop();CHKERRQ(ierr); 2055 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error"); 2056 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */ 2057 2058 /* Extract an explicit representation of Q */ 2059 Q = Ainv; 2060 ierr = PetscArraycpy(Q,A,mstride*n);CHKERRQ(ierr); 2061 K = N; /* full rank */ 2062 PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info)); 2063 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error"); 2064 2065 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */ 2066 Alpha = 1.0; 2067 ldb = lda; 2068 PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb)); 2069 /* Ainv is Q, overwritten with inverse */ 2070 2071 #if defined(PETSC_USE_COMPLEX) 2072 { 2073 PetscInt i; 2074 for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]); 2075 ierr = PetscFree2(A,Ainv);CHKERRQ(ierr); 2076 } 2077 #endif 2078 PetscFunctionReturn(0); 2079 } 2080 2081 /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */ 2082 static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B) 2083 { 2084 PetscErrorCode ierr; 2085 PetscReal *Bv; 2086 PetscInt i,j; 2087 2088 PetscFunctionBegin; 2089 ierr = PetscMalloc1((ninterval+1)*ndegree,&Bv);CHKERRQ(ierr); 2090 /* Point evaluation of L_p on all the source vertices */ 2091 ierr = PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);CHKERRQ(ierr); 2092 /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */ 2093 for (i=0; i<ninterval; i++) { 2094 for (j=0; j<ndegree; j++) { 2095 if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j]; 2096 else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j]; 2097 } 2098 } 2099 ierr = PetscFree(Bv);CHKERRQ(ierr); 2100 PetscFunctionReturn(0); 2101 } 2102 2103 /*@ 2104 PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals 2105 2106 Not Collective 2107 2108 Input Parameters: 2109 + degree - degree of reconstruction polynomial 2110 . nsource - number of source intervals 2111 . sourcex - sorted coordinates of source cell boundaries (length nsource+1) 2112 . ntarget - number of target intervals 2113 - targetx - sorted coordinates of target cell boundaries (length ntarget+1) 2114 2115 Output Parameter: 2116 . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s] 2117 2118 Level: advanced 2119 2120 .seealso: PetscDTLegendreEval() 2121 @*/ 2122 PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R) 2123 { 2124 PetscErrorCode ierr; 2125 PetscInt i,j,k,*bdegrees,worksize; 2126 PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget; 2127 PetscScalar *tau,*work; 2128 2129 PetscFunctionBegin; 2130 PetscValidRealPointer(sourcex,3); 2131 PetscValidRealPointer(targetx,5); 2132 PetscValidRealPointer(R,6); 2133 if (degree >= nsource) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Reconstruction degree %D must be less than number of source intervals %D",degree,nsource); 2134 if (PetscDefined(USE_DEBUG)) { 2135 for (i=0; i<nsource; i++) { 2136 if (sourcex[i] >= sourcex[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Source interval %D has negative orientation (%g,%g)",i,(double)sourcex[i],(double)sourcex[i+1]); 2137 } 2138 for (i=0; i<ntarget; i++) { 2139 if (targetx[i] >= targetx[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Target interval %D has negative orientation (%g,%g)",i,(double)targetx[i],(double)targetx[i+1]); 2140 } 2141 } 2142 xmin = PetscMin(sourcex[0],targetx[0]); 2143 xmax = PetscMax(sourcex[nsource],targetx[ntarget]); 2144 center = (xmin + xmax)/2; 2145 hscale = (xmax - xmin)/2; 2146 worksize = nsource; 2147 ierr = PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);CHKERRQ(ierr); 2148 ierr = PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);CHKERRQ(ierr); 2149 for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale; 2150 for (i=0; i<=degree; i++) bdegrees[i] = i+1; 2151 ierr = PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);CHKERRQ(ierr); 2152 ierr = PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);CHKERRQ(ierr); 2153 for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale; 2154 ierr = PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);CHKERRQ(ierr); 2155 for (i=0; i<ntarget; i++) { 2156 PetscReal rowsum = 0; 2157 for (j=0; j<nsource; j++) { 2158 PetscReal sum = 0; 2159 for (k=0; k<degree+1; k++) { 2160 sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j]; 2161 } 2162 R[i*nsource+j] = sum; 2163 rowsum += sum; 2164 } 2165 for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */ 2166 } 2167 ierr = PetscFree4(bdegrees,sourcey,Bsource,work);CHKERRQ(ierr); 2168 ierr = PetscFree4(tau,Bsinv,targety,Btarget);CHKERRQ(ierr); 2169 PetscFunctionReturn(0); 2170 } 2171 2172 /*@C 2173 PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points 2174 2175 Not Collective 2176 2177 Input Parameters: 2178 + n - the number of GLL nodes 2179 . nodes - the GLL nodes 2180 . weights - the GLL weights 2181 - f - the function values at the nodes 2182 2183 Output Parameter: 2184 . in - the value of the integral 2185 2186 Level: beginner 2187 2188 .seealso: PetscDTGaussLobattoLegendreQuadrature() 2189 2190 @*/ 2191 PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal *nodes,PetscReal *weights,const PetscReal *f,PetscReal *in) 2192 { 2193 PetscInt i; 2194 2195 PetscFunctionBegin; 2196 *in = 0.; 2197 for (i=0; i<n; i++) { 2198 *in += f[i]*f[i]*weights[i]; 2199 } 2200 PetscFunctionReturn(0); 2201 } 2202 2203 /*@C 2204 PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element 2205 2206 Not Collective 2207 2208 Input Parameters: 2209 + n - the number of GLL nodes 2210 . nodes - the GLL nodes 2211 - weights - the GLL weights 2212 2213 Output Parameter: 2214 . A - the stiffness element 2215 2216 Level: beginner 2217 2218 Notes: 2219 Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy() 2220 2221 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) 2222 2223 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy() 2224 2225 @*/ 2226 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 2227 { 2228 PetscReal **A; 2229 PetscErrorCode ierr; 2230 const PetscReal *gllnodes = nodes; 2231 const PetscInt p = n-1; 2232 PetscReal z0,z1,z2 = -1,x,Lpj,Lpr; 2233 PetscInt i,j,nn,r; 2234 2235 PetscFunctionBegin; 2236 ierr = PetscMalloc1(n,&A);CHKERRQ(ierr); 2237 ierr = PetscMalloc1(n*n,&A[0]);CHKERRQ(ierr); 2238 for (i=1; i<n; i++) A[i] = A[i-1]+n; 2239 2240 for (j=1; j<p; j++) { 2241 x = gllnodes[j]; 2242 z0 = 1.; 2243 z1 = x; 2244 for (nn=1; nn<p; nn++) { 2245 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 2246 z0 = z1; 2247 z1 = z2; 2248 } 2249 Lpj=z2; 2250 for (r=1; r<p; r++) { 2251 if (r == j) { 2252 A[j][j]=2./(3.*(1.-gllnodes[j]*gllnodes[j])*Lpj*Lpj); 2253 } else { 2254 x = gllnodes[r]; 2255 z0 = 1.; 2256 z1 = x; 2257 for (nn=1; nn<p; nn++) { 2258 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 2259 z0 = z1; 2260 z1 = z2; 2261 } 2262 Lpr = z2; 2263 A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(gllnodes[j]-gllnodes[r])*(gllnodes[j]-gllnodes[r])); 2264 } 2265 } 2266 } 2267 for (j=1; j<p+1; j++) { 2268 x = gllnodes[j]; 2269 z0 = 1.; 2270 z1 = x; 2271 for (nn=1; nn<p; nn++) { 2272 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 2273 z0 = z1; 2274 z1 = z2; 2275 } 2276 Lpj = z2; 2277 A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+gllnodes[j])*(1.+gllnodes[j])); 2278 A[0][j] = A[j][0]; 2279 } 2280 for (j=0; j<p; j++) { 2281 x = gllnodes[j]; 2282 z0 = 1.; 2283 z1 = x; 2284 for (nn=1; nn<p; nn++) { 2285 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 2286 z0 = z1; 2287 z1 = z2; 2288 } 2289 Lpj=z2; 2290 2291 A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-gllnodes[j])*(1.-gllnodes[j])); 2292 A[j][p] = A[p][j]; 2293 } 2294 A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.; 2295 A[p][p]=A[0][0]; 2296 *AA = A; 2297 PetscFunctionReturn(0); 2298 } 2299 2300 /*@C 2301 PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element 2302 2303 Not Collective 2304 2305 Input Parameters: 2306 + n - the number of GLL nodes 2307 . nodes - the GLL nodes 2308 . weights - the GLL weightss 2309 - A - the stiffness element 2310 2311 Level: beginner 2312 2313 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate() 2314 2315 @*/ 2316 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 2317 { 2318 PetscErrorCode ierr; 2319 2320 PetscFunctionBegin; 2321 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 2322 ierr = PetscFree(*AA);CHKERRQ(ierr); 2323 *AA = NULL; 2324 PetscFunctionReturn(0); 2325 } 2326 2327 /*@C 2328 PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element 2329 2330 Not Collective 2331 2332 Input Parameter: 2333 + n - the number of GLL nodes 2334 . nodes - the GLL nodes 2335 . weights - the GLL weights 2336 2337 Output Parameters: 2338 . AA - the stiffness element 2339 - AAT - the transpose of AA (pass in NULL if you do not need this array) 2340 2341 Level: beginner 2342 2343 Notes: 2344 Destroy this with PetscGaussLobattoLegendreElementGradientDestroy() 2345 2346 You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented 2347 2348 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy() 2349 2350 @*/ 2351 PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT) 2352 { 2353 PetscReal **A, **AT = NULL; 2354 PetscErrorCode ierr; 2355 const PetscReal *gllnodes = nodes; 2356 const PetscInt p = n-1; 2357 PetscReal Li, Lj,d0; 2358 PetscInt i,j; 2359 2360 PetscFunctionBegin; 2361 ierr = PetscMalloc1(n,&A);CHKERRQ(ierr); 2362 ierr = PetscMalloc1(n*n,&A[0]);CHKERRQ(ierr); 2363 for (i=1; i<n; i++) A[i] = A[i-1]+n; 2364 2365 if (AAT) { 2366 ierr = PetscMalloc1(n,&AT);CHKERRQ(ierr); 2367 ierr = PetscMalloc1(n*n,&AT[0]);CHKERRQ(ierr); 2368 for (i=1; i<n; i++) AT[i] = AT[i-1]+n; 2369 } 2370 2371 if (n==1) {A[0][0] = 0.;} 2372 d0 = (PetscReal)p*((PetscReal)p+1.)/4.; 2373 for (i=0; i<n; i++) { 2374 for (j=0; j<n; j++) { 2375 A[i][j] = 0.; 2376 ierr = PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li);CHKERRQ(ierr); 2377 ierr = PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj);CHKERRQ(ierr); 2378 if (i!=j) A[i][j] = Li/(Lj*(gllnodes[i]-gllnodes[j])); 2379 if ((j==i) && (i==0)) A[i][j] = -d0; 2380 if (j==i && i==p) A[i][j] = d0; 2381 if (AT) AT[j][i] = A[i][j]; 2382 } 2383 } 2384 if (AAT) *AAT = AT; 2385 *AA = A; 2386 PetscFunctionReturn(0); 2387 } 2388 2389 /*@C 2390 PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGaussLobattoLegendreElementGradientCreate() 2391 2392 Not Collective 2393 2394 Input Parameters: 2395 + n - the number of GLL nodes 2396 . nodes - the GLL nodes 2397 . weights - the GLL weights 2398 . AA - the stiffness element 2399 - AAT - the transpose of the element 2400 2401 Level: beginner 2402 2403 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionCreate() 2404 2405 @*/ 2406 PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT) 2407 { 2408 PetscErrorCode ierr; 2409 2410 PetscFunctionBegin; 2411 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 2412 ierr = PetscFree(*AA);CHKERRQ(ierr); 2413 *AA = NULL; 2414 if (*AAT) { 2415 ierr = PetscFree((*AAT)[0]);CHKERRQ(ierr); 2416 ierr = PetscFree(*AAT);CHKERRQ(ierr); 2417 *AAT = NULL; 2418 } 2419 PetscFunctionReturn(0); 2420 } 2421 2422 /*@C 2423 PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element 2424 2425 Not Collective 2426 2427 Input Parameters: 2428 + n - the number of GLL nodes 2429 . nodes - the GLL nodes 2430 - weights - the GLL weightss 2431 2432 Output Parameter: 2433 . AA - the stiffness element 2434 2435 Level: beginner 2436 2437 Notes: 2438 Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy() 2439 2440 This is the same as the Gradient operator multiplied by the diagonal mass matrix 2441 2442 You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented 2443 2444 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionDestroy() 2445 2446 @*/ 2447 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 2448 { 2449 PetscReal **D; 2450 PetscErrorCode ierr; 2451 const PetscReal *gllweights = weights; 2452 const PetscInt glln = n; 2453 PetscInt i,j; 2454 2455 PetscFunctionBegin; 2456 ierr = PetscGaussLobattoLegendreElementGradientCreate(n,nodes,weights,&D,NULL);CHKERRQ(ierr); 2457 for (i=0; i<glln; i++) { 2458 for (j=0; j<glln; j++) { 2459 D[i][j] = gllweights[i]*D[i][j]; 2460 } 2461 } 2462 *AA = D; 2463 PetscFunctionReturn(0); 2464 } 2465 2466 /*@C 2467 PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element 2468 2469 Not Collective 2470 2471 Input Parameters: 2472 + n - the number of GLL nodes 2473 . nodes - the GLL nodes 2474 . weights - the GLL weights 2475 - A - advection 2476 2477 Level: beginner 2478 2479 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementAdvectionCreate() 2480 2481 @*/ 2482 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 2483 { 2484 PetscErrorCode ierr; 2485 2486 PetscFunctionBegin; 2487 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 2488 ierr = PetscFree(*AA);CHKERRQ(ierr); 2489 *AA = NULL; 2490 PetscFunctionReturn(0); 2491 } 2492 2493 PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 2494 { 2495 PetscReal **A; 2496 PetscErrorCode ierr; 2497 const PetscReal *gllweights = weights; 2498 const PetscInt glln = n; 2499 PetscInt i,j; 2500 2501 PetscFunctionBegin; 2502 ierr = PetscMalloc1(glln,&A);CHKERRQ(ierr); 2503 ierr = PetscMalloc1(glln*glln,&A[0]);CHKERRQ(ierr); 2504 for (i=1; i<glln; i++) A[i] = A[i-1]+glln; 2505 if (glln==1) {A[0][0] = 0.;} 2506 for (i=0; i<glln; i++) { 2507 for (j=0; j<glln; j++) { 2508 A[i][j] = 0.; 2509 if (j==i) A[i][j] = gllweights[i]; 2510 } 2511 } 2512 *AA = A; 2513 PetscFunctionReturn(0); 2514 } 2515 2516 PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 2517 { 2518 PetscErrorCode ierr; 2519 2520 PetscFunctionBegin; 2521 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 2522 ierr = PetscFree(*AA);CHKERRQ(ierr); 2523 *AA = NULL; 2524 PetscFunctionReturn(0); 2525 } 2526 2527 /*@ 2528 PetscDTIndexToBary - convert an index into a barycentric coordinate. 2529 2530 Input Parameters: 2531 + 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) 2532 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 2533 - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum) 2534 2535 Output Parameter: 2536 . coord - will be filled with the barycentric coordinate 2537 2538 Level: beginner 2539 2540 Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the 2541 least significant and the last index is the most significant. 2542 2543 .seealso: PetscDTBaryToIndex() 2544 @*/ 2545 PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[]) 2546 { 2547 PetscInt c, d, s, total, subtotal, nexttotal; 2548 2549 PetscFunctionBeginHot; 2550 if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 2551 if (index < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative"); 2552 if (!len) { 2553 if (!sum && !index) PetscFunctionReturn(0); 2554 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 2555 } 2556 for (c = 1, total = 1; c <= len; c++) { 2557 /* total is the number of ways to have a tuple of length c with sum */ 2558 if (index < total) break; 2559 total = (total * (sum + c)) / c; 2560 } 2561 if (c > len) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range"); 2562 for (d = c; d < len; d++) coord[d] = 0; 2563 for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) { 2564 /* subtotal is the number of ways to have a tuple of length c with sum s */ 2565 /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */ 2566 if ((index + subtotal) >= total) { 2567 coord[--c] = sum - s; 2568 index -= (total - subtotal); 2569 sum = s; 2570 total = nexttotal; 2571 subtotal = 1; 2572 nexttotal = 1; 2573 s = 0; 2574 } else { 2575 subtotal = (subtotal * (c + s)) / (s + 1); 2576 nexttotal = (nexttotal * (c - 1 + s)) / (s + 1); 2577 s++; 2578 } 2579 } 2580 PetscFunctionReturn(0); 2581 } 2582 2583 /*@ 2584 PetscDTBaryToIndex - convert a barycentric coordinate to an index 2585 2586 Input Parameters: 2587 + 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) 2588 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 2589 - coord - a barycentric coordinate with the given length and sum 2590 2591 Output Parameter: 2592 . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum) 2593 2594 Level: beginner 2595 2596 Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the 2597 least significant and the last index is the most significant. 2598 2599 .seealso: PetscDTIndexToBary 2600 @*/ 2601 PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index) 2602 { 2603 PetscInt c; 2604 PetscInt i; 2605 PetscInt total; 2606 2607 PetscFunctionBeginHot; 2608 if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 2609 if (!len) { 2610 if (!sum) { 2611 *index = 0; 2612 PetscFunctionReturn(0); 2613 } 2614 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 2615 } 2616 for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c; 2617 i = total - 1; 2618 c = len - 1; 2619 sum -= coord[c]; 2620 while (sum > 0) { 2621 PetscInt subtotal; 2622 PetscInt s; 2623 2624 for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s; 2625 i -= subtotal; 2626 sum -= coord[--c]; 2627 } 2628 *index = i; 2629 PetscFunctionReturn(0); 2630 } 2631