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_", 0}; 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, 2); 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 Arguments: 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 Arguments: 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, 4); 480 q->points = points; 481 } 482 if (weights) { 483 PetscValidPointer(weights, 5); 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 Parameter: 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 static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p) 632 { 633 PetscReal ak, bk; 634 PetscReal abk1; 635 PetscInt i,l,maxdegree; 636 637 PetscFunctionBegin; 638 maxdegree = degrees[ndegree-1] - k; 639 ak = a + k; 640 bk = b + k; 641 abk1 = a + b + k + 1.; 642 if (maxdegree < 0) { 643 for (i = 0; i < npoints; i++) for (l = 0; l < ndegree; l++) p[i*ndegree+l] = 0.; 644 PetscFunctionReturn(0); 645 } 646 for (i=0; i<npoints; i++) { 647 PetscReal pm1,pm2,x; 648 PetscReal cnm1, cnm1x, cnm2; 649 PetscInt j,m; 650 651 x = points[i]; 652 pm2 = 1.; 653 PetscDTJacobiRecurrence_Internal(1,ak,bk,cnm1,cnm1x,cnm2); 654 pm1 = (cnm1 + cnm1x*x); 655 l = 0; 656 while (l < ndegree && degrees[l] - k < 0) { 657 p[l++] = 0.; 658 } 659 while (l < ndegree && degrees[l] - k == 0) { 660 p[l] = pm2; 661 for (m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5; 662 l++; 663 } 664 while (l < ndegree && degrees[l] - k == 1) { 665 p[l] = pm1; 666 for (m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5; 667 l++; 668 } 669 for (j=2; j<=maxdegree; j++) { 670 PetscReal pp; 671 672 PetscDTJacobiRecurrence_Internal(j,ak,bk,cnm1,cnm1x,cnm2); 673 pp = (cnm1 + cnm1x*x)*pm1 - cnm2*pm2; 674 pm2 = pm1; 675 pm1 = pp; 676 while (l < ndegree && degrees[l] - k == j) { 677 p[l] = pp; 678 for (m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5; 679 l++; 680 } 681 } 682 p += ndegree; 683 } 684 PetscFunctionReturn(0); 685 } 686 687 /*@ 688 PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$ 689 at points 690 691 Not Collective 692 693 Input Arguments: 694 + npoints - number of spatial points to evaluate at 695 . alpha - the left exponent > -1 696 . beta - the right exponent > -1 697 . points - array of locations to evaluate at 698 . ndegree - number of basis degrees to evaluate 699 - degrees - sorted array of degrees to evaluate 700 701 Output Arguments: 702 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL) 703 . D - row-oriented derivative evaluation matrix (or NULL) 704 - D2 - row-oriented second derivative evaluation matrix (or NULL) 705 706 Level: intermediate 707 708 .seealso: PetscDTGaussQuadrature() 709 @*/ 710 PetscErrorCode PetscDTJacobiEval(PetscInt npoints,PetscReal alpha, PetscReal beta, const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2) 711 { 712 PetscErrorCode ierr; 713 714 PetscFunctionBegin; 715 if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1."); 716 if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1."); 717 if (!npoints || !ndegree) PetscFunctionReturn(0); 718 if (B) {ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B);CHKERRQ(ierr);} 719 if (D) {ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D);CHKERRQ(ierr);} 720 if (D2) {ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2);CHKERRQ(ierr);} 721 PetscFunctionReturn(0); 722 } 723 724 /*@ 725 PetscDTLegendreEval - evaluate Legendre polynomials at points 726 727 Not Collective 728 729 Input Arguments: 730 + npoints - number of spatial points to evaluate at 731 . points - array of locations to evaluate at 732 . ndegree - number of basis degrees to evaluate 733 - degrees - sorted array of degrees to evaluate 734 735 Output Arguments: 736 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL) 737 . D - row-oriented derivative evaluation matrix (or NULL) 738 - D2 - row-oriented second derivative evaluation matrix (or NULL) 739 740 Level: intermediate 741 742 .seealso: PetscDTGaussQuadrature() 743 @*/ 744 PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2) 745 { 746 PetscErrorCode ierr; 747 748 PetscFunctionBegin; 749 ierr = PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2);CHKERRQ(ierr); 750 PetscFunctionReturn(0); 751 } 752 753 /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V 754 * with lds n; diag and subdiag are overwritten */ 755 static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[], 756 PetscReal eigs[], PetscScalar V[]) 757 { 758 char jobz = 'V'; /* eigenvalues and eigenvectors */ 759 char range = 'A'; /* all eigenvalues will be found */ 760 PetscReal VL = 0.; /* ignored because range is 'A' */ 761 PetscReal VU = 0.; /* ignored because range is 'A' */ 762 PetscBLASInt IL = 0; /* ignored because range is 'A' */ 763 PetscBLASInt IU = 0; /* ignored because range is 'A' */ 764 PetscReal abstol = 0.; /* unused */ 765 PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */ 766 PetscBLASInt *isuppz; 767 PetscBLASInt lwork, liwork; 768 PetscReal workquery; 769 PetscBLASInt iworkquery; 770 PetscBLASInt *iwork; 771 PetscBLASInt info; 772 PetscReal *work = NULL; 773 PetscErrorCode ierr; 774 775 PetscFunctionBegin; 776 #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG) 777 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found"); 778 #endif 779 ierr = PetscBLASIntCast(n, &bn);CHKERRQ(ierr); 780 ierr = PetscBLASIntCast(n, &ldz);CHKERRQ(ierr); 781 #if !defined(PETSC_MISSING_LAPACK_STEGR) 782 ierr = PetscMalloc1(2 * n, &isuppz);CHKERRQ(ierr); 783 lwork = -1; 784 liwork = -1; 785 PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,&workquery,&lwork,&iworkquery,&liwork,&info)); 786 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEGR error"); 787 lwork = (PetscBLASInt) workquery; 788 liwork = (PetscBLASInt) iworkquery; 789 ierr = PetscMalloc2(lwork, &work, liwork, &iwork);CHKERRQ(ierr); 790 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 791 PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,work,&lwork,iwork,&liwork,&info)); 792 ierr = PetscFPTrapPop();CHKERRQ(ierr); 793 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEGR error"); 794 ierr = PetscFree2(work, iwork);CHKERRQ(ierr); 795 ierr = PetscFree(isuppz);CHKERRQ(ierr); 796 #elif !defined(PETSC_MISSING_LAPACK_STEQR) 797 jobz = 'I'; /* Compute eigenvalues and eigenvectors of the 798 tridiagonal matrix. Z is initialized to the identity 799 matrix. */ 800 ierr = PetscMalloc1(PetscMax(1,2*n-2),&work);CHKERRQ(ierr); 801 PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&bn,diag,subdiag,V,&ldz,work,&info)); 802 ierr = PetscFPTrapPop();CHKERRQ(ierr); 803 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error"); 804 ierr = PetscFree(work);CHKERRQ(ierr); 805 ierr = PetscArraycpy(eigs,diag,n);CHKERRQ(ierr); 806 #endif 807 PetscFunctionReturn(0); 808 } 809 810 /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi 811 * quadrature rules on the interval [-1, 1] */ 812 static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw) 813 { 814 PetscReal twoab1; 815 PetscInt m = n - 2; 816 PetscReal a = alpha + 1.; 817 PetscReal b = beta + 1.; 818 PetscReal gra, grb; 819 820 PetscFunctionBegin; 821 twoab1 = PetscPowReal(2., a + b - 1.); 822 #if defined(PETSC_HAVE_LGAMMA) 823 grb = PetscExpReal(2. * PetscLGamma(b+1.) + PetscLGamma(m+1.) + PetscLGamma(m+a+1.) - 824 (PetscLGamma(m+b+1) + PetscLGamma(m+a+b+1.))); 825 gra = PetscExpReal(2. * PetscLGamma(a+1.) + PetscLGamma(m+1.) + PetscLGamma(m+b+1.) - 826 (PetscLGamma(m+a+1) + PetscLGamma(m+a+b+1.))); 827 #else 828 { 829 PetscInt alphai = (PetscInt) alpha; 830 PetscInt betai = (PetscInt) beta; 831 PetscErrorCode ierr; 832 833 if ((PetscReal) alphai == alpha && (PetscReal) betai == beta) { 834 PetscReal binom1, binom2; 835 836 ierr = PetscDTBinomial(m+b, b, &binom1);CHKERRQ(ierr); 837 ierr = PetscDTBinomial(m+a+b, b, &binom2);CHKERRQ(ierr); 838 grb = 1./ (binom1 * binom2); 839 ierr = PetscDTBinomial(m+a, a, &binom1);CHKERRQ(ierr); 840 ierr = PetscDTBinomial(m+a+b, a, &binom2);CHKERRQ(ierr); 841 gra = 1./ (binom1 * binom2); 842 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable."); 843 } 844 #endif 845 *leftw = twoab1 * grb / b; 846 *rightw = twoab1 * gra / a; 847 PetscFunctionReturn(0); 848 } 849 850 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x. 851 Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */ 852 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 853 { 854 PetscReal pn1, pn2; 855 PetscReal cnm1, cnm1x, cnm2; 856 PetscInt k; 857 858 PetscFunctionBegin; 859 if (!n) {*P = 1.0; PetscFunctionReturn(0);} 860 PetscDTJacobiRecurrence_Internal(1,a,b,cnm1,cnm1x,cnm2); 861 pn2 = 1.; 862 pn1 = cnm1 + cnm1x*x; 863 if (n == 1) {*P = pn1; PetscFunctionReturn(0);} 864 *P = 0.0; 865 for (k = 2; k < n+1; ++k) { 866 PetscDTJacobiRecurrence_Internal(k,a,b,cnm1,cnm1x,cnm2); 867 868 *P = (cnm1 + cnm1x*x)*pn1 - cnm2*pn2; 869 pn2 = pn1; 870 pn1 = *P; 871 } 872 PetscFunctionReturn(0); 873 } 874 875 /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */ 876 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P) 877 { 878 PetscReal nP; 879 PetscInt i; 880 PetscErrorCode ierr; 881 882 PetscFunctionBegin; 883 *P = 0.0; 884 if (k > n) PetscFunctionReturn(0); 885 ierr = PetscDTComputeJacobi(a+k, b+k, n-k, x, &nP);CHKERRQ(ierr); 886 for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5; 887 *P = nP; 888 PetscFunctionReturn(0); 889 } 890 891 /* Maps from [-1,1]^2 to the (-1,1) reference triangle */ 892 PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta) 893 { 894 PetscFunctionBegin; 895 *xi = 0.5 * (1.0 + x) * (1.0 - y) - 1.0; 896 *eta = y; 897 PetscFunctionReturn(0); 898 } 899 900 static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[]) 901 { 902 PetscInt maxIter = 100; 903 PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON)); 904 PetscReal a1, a6, gf; 905 PetscInt k; 906 PetscErrorCode ierr; 907 908 PetscFunctionBegin; 909 910 a1 = PetscPowReal(2.0, a+b+1); 911 #if defined(PETSC_HAVE_LGAMMA) 912 { 913 PetscReal a2, a3, a4, a5; 914 a2 = PetscLGamma(a + npoints + 1); 915 a3 = PetscLGamma(b + npoints + 1); 916 a4 = PetscLGamma(a + b + npoints + 1); 917 a5 = PetscLGamma(npoints + 1); 918 gf = PetscExpReal(a2 + a3 - (a4 + a5)); 919 } 920 #else 921 { 922 PetscInt ia, ib; 923 924 ia = (PetscInt) a; 925 ib = (PetscInt) b; 926 gf = 1.; 927 if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */ 928 for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k); 929 } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */ 930 for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k); 931 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable."); 932 } 933 #endif 934 935 a6 = a1 * gf; 936 /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses. 937 Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */ 938 for (k = 0; k < npoints; ++k) { 939 PetscReal r = PetscCosReal(PETSC_PI * (1. - (4.*k + 3. + 2.*b) / (4.*npoints + 2.*(a + b + 1.)))), dP; 940 PetscInt j; 941 942 if (k > 0) r = 0.5 * (r + x[k-1]); 943 for (j = 0; j < maxIter; ++j) { 944 PetscReal s = 0.0, delta, f, fp; 945 PetscInt i; 946 947 for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]); 948 ierr = PetscDTComputeJacobi(a, b, npoints, r, &f);CHKERRQ(ierr); 949 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp);CHKERRQ(ierr); 950 delta = f / (fp - f * s); 951 r = r - delta; 952 if (PetscAbsReal(delta) < eps) break; 953 } 954 x[k] = r; 955 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP);CHKERRQ(ierr); 956 w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP); 957 } 958 PetscFunctionReturn(0); 959 } 960 961 /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi 962 * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */ 963 static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s) 964 { 965 PetscInt i; 966 967 PetscFunctionBegin; 968 for (i = 0; i < nPoints; i++) { 969 PetscReal A, B, C; 970 971 PetscDTJacobiRecurrence_Internal(i+1,a,b,A,B,C); 972 d[i] = -A / B; 973 if (i) s[i-1] *= C / B; 974 if (i < nPoints - 1) s[i] = 1. / B; 975 } 976 PetscFunctionReturn(0); 977 } 978 979 static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w) 980 { 981 PetscReal mu0; 982 PetscReal ga, gb, gab; 983 PetscInt i; 984 PetscErrorCode ierr; 985 986 PetscFunctionBegin; 987 ierr = PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite);CHKERRQ(ierr); 988 989 #if defined(PETSC_HAVE_TGAMMA) 990 ga = PetscTGamma(a + 1); 991 gb = PetscTGamma(b + 1); 992 gab = PetscTGamma(a + b + 2); 993 #else 994 { 995 PetscInt ia, ib; 996 997 ia = (PetscInt) a; 998 ib = (PetscInt) b; 999 if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */ 1000 ierr = PetscDTFactorial(ia, &ga);CHKERRQ(ierr); 1001 ierr = PetscDTFactorial(ib, &gb);CHKERRQ(ierr); 1002 ierr = PetscDTFactorial(ia + ib + 1, &gb);CHKERRQ(ierr); 1003 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable."); 1004 } 1005 #endif 1006 mu0 = PetscPowReal(2.,a + b + 1.) * ga * gb / gab; 1007 1008 #if defined(PETSCDTGAUSSIANQUADRATURE_EIG) 1009 { 1010 PetscReal *diag, *subdiag; 1011 PetscScalar *V; 1012 1013 ierr = PetscMalloc2(npoints, &diag, npoints, &subdiag);CHKERRQ(ierr); 1014 ierr = PetscMalloc1(npoints*npoints, &V);CHKERRQ(ierr); 1015 ierr = PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag);CHKERRQ(ierr); 1016 for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]); 1017 ierr = PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V);CHKERRQ(ierr); 1018 for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0; 1019 ierr = PetscFree(V);CHKERRQ(ierr); 1020 ierr = PetscFree2(diag, subdiag);CHKERRQ(ierr); 1021 } 1022 #else 1023 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found"); 1024 #endif 1025 { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the 1026 eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that 1027 the eigenvalues are sorted */ 1028 PetscBool sorted; 1029 1030 ierr = PetscSortedReal(npoints, x, &sorted);CHKERRQ(ierr); 1031 if (!sorted) { 1032 PetscInt *order, i; 1033 PetscReal *tmp; 1034 1035 ierr = PetscMalloc2(npoints, &order, npoints, &tmp);CHKERRQ(ierr); 1036 for (i = 0; i < npoints; i++) order[i] = i; 1037 ierr = PetscSortRealWithPermutation(npoints, x, order);CHKERRQ(ierr); 1038 ierr = PetscArraycpy(tmp, x, npoints);CHKERRQ(ierr); 1039 for (i = 0; i < npoints; i++) x[i] = tmp[order[i]]; 1040 ierr = PetscArraycpy(tmp, w, npoints);CHKERRQ(ierr); 1041 for (i = 0; i < npoints; i++) w[i] = tmp[order[i]]; 1042 ierr = PetscFree2(order, tmp);CHKERRQ(ierr); 1043 } 1044 } 1045 PetscFunctionReturn(0); 1046 } 1047 1048 static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1049 { 1050 PetscErrorCode ierr; 1051 1052 PetscFunctionBegin; 1053 if (npoints < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be positive"); 1054 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1055 if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1."); 1056 if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1."); 1057 1058 if (newton) { 1059 ierr = PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w);CHKERRQ(ierr); 1060 } else { 1061 ierr = PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w);CHKERRQ(ierr); 1062 } 1063 if (alpha == beta) { /* symmetrize */ 1064 PetscInt i; 1065 for (i = 0; i < (npoints + 1) / 2; i++) { 1066 PetscInt j = npoints - 1 - i; 1067 PetscReal xi = x[i]; 1068 PetscReal xj = x[j]; 1069 PetscReal wi = w[i]; 1070 PetscReal wj = w[j]; 1071 1072 x[i] = (xi - xj) / 2.; 1073 x[j] = (xj - xi) / 2.; 1074 w[i] = w[j] = (wi + wj) / 2.; 1075 } 1076 } 1077 PetscFunctionReturn(0); 1078 } 1079 1080 /*@ 1081 PetscDTGaussJacobiQuadrature - quadrature for the interval [a, b] with the weight function 1082 $(x-a)^\alpha (x-b)^\beta$. 1083 1084 Not collective 1085 1086 Input Parameters: 1087 + npoints - the number of points in the quadrature rule 1088 . a - the left endpoint of the interval 1089 . b - the right endpoint of the interval 1090 . alpha - the left exponent 1091 - beta - the right exponent 1092 1093 Output Parameters: 1094 + x - array of length npoints, the locations of the quadrature points 1095 - w - array of length npoints, the weights of the quadrature points 1096 1097 Level: intermediate 1098 1099 Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 1. 1100 @*/ 1101 PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1102 { 1103 PetscInt i; 1104 PetscErrorCode ierr; 1105 1106 PetscFunctionBegin; 1107 ierr = PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);CHKERRQ(ierr); 1108 if (a != -1. || b != 1.) { /* shift */ 1109 for (i = 0; i < npoints; i++) { 1110 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1111 w[i] *= (b - a) / 2.; 1112 } 1113 } 1114 PetscFunctionReturn(0); 1115 } 1116 1117 static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton) 1118 { 1119 PetscInt i; 1120 PetscErrorCode ierr; 1121 1122 PetscFunctionBegin; 1123 if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be positive"); 1124 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */ 1125 if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1."); 1126 if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1."); 1127 1128 x[0] = -1.; 1129 x[npoints-1] = 1.; 1130 if (npoints > 2) { 1131 ierr = PetscDTGaussJacobiQuadrature_Internal(npoints-2, alpha+1., beta+1., &x[1], &w[1], newton);CHKERRQ(ierr); 1132 } 1133 for (i = 1; i < npoints - 1; i++) { 1134 w[i] /= (1. - x[i]*x[i]); 1135 } 1136 ierr = PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints-1]);CHKERRQ(ierr); 1137 PetscFunctionReturn(0); 1138 } 1139 1140 /*@ 1141 PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval [a, b] with the weight function 1142 $(x-a)^\alpha (x-b)^\beta$, with endpoints a and b included as quadrature points. 1143 1144 Not collective 1145 1146 Input Parameters: 1147 + npoints - the number of points in the quadrature rule 1148 . a - the left endpoint of the interval 1149 . b - the right endpoint of the interval 1150 . alpha - the left exponent 1151 - beta - the right exponent 1152 1153 Output Parameters: 1154 + x - array of length npoints, the locations of the quadrature points 1155 - w - array of length npoints, the weights of the quadrature points 1156 1157 Level: intermediate 1158 1159 Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 3. 1160 @*/ 1161 PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[]) 1162 { 1163 PetscInt i; 1164 PetscErrorCode ierr; 1165 1166 PetscFunctionBegin; 1167 ierr = PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);CHKERRQ(ierr); 1168 if (a != -1. || b != 1.) { /* shift */ 1169 for (i = 0; i < npoints; i++) { 1170 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1171 w[i] *= (b - a) / 2.; 1172 } 1173 } 1174 PetscFunctionReturn(0); 1175 } 1176 1177 /*@ 1178 PetscDTGaussQuadrature - create Gauss-Legendre quadrature 1179 1180 Not Collective 1181 1182 Input Arguments: 1183 + npoints - number of points 1184 . a - left end of interval (often-1) 1185 - b - right end of interval (often +1) 1186 1187 Output Arguments: 1188 + x - quadrature points 1189 - w - quadrature weights 1190 1191 Level: intermediate 1192 1193 References: 1194 . 1. - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969. 1195 1196 .seealso: PetscDTLegendreEval() 1197 @*/ 1198 PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w) 1199 { 1200 PetscInt i; 1201 PetscErrorCode ierr; 1202 1203 PetscFunctionBegin; 1204 ierr = PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal);CHKERRQ(ierr); 1205 if (a != -1. || b != 1.) { /* shift */ 1206 for (i = 0; i < npoints; i++) { 1207 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a; 1208 w[i] *= (b - a) / 2.; 1209 } 1210 } 1211 PetscFunctionReturn(0); 1212 } 1213 1214 /*@C 1215 PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre 1216 nodes of a given size on the domain [-1,1] 1217 1218 Not Collective 1219 1220 Input Parameter: 1221 + n - number of grid nodes 1222 - type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON 1223 1224 Output Arguments: 1225 + x - quadrature points 1226 - w - quadrature weights 1227 1228 Notes: 1229 For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not 1230 close enough to the desired solution 1231 1232 These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes 1233 1234 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 1235 1236 Level: intermediate 1237 1238 .seealso: PetscDTGaussQuadrature() 1239 1240 @*/ 1241 PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal *x,PetscReal *w) 1242 { 1243 PetscBool newton; 1244 PetscErrorCode ierr; 1245 1246 PetscFunctionBegin; 1247 if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide at least 2 grid points per element"); 1248 newton = (PetscBool) (type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON); 1249 ierr = PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton);CHKERRQ(ierr); 1250 PetscFunctionReturn(0); 1251 } 1252 1253 /*@ 1254 PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature 1255 1256 Not Collective 1257 1258 Input Arguments: 1259 + dim - The spatial dimension 1260 . Nc - The number of components 1261 . npoints - number of points in one dimension 1262 . a - left end of interval (often-1) 1263 - b - right end of interval (often +1) 1264 1265 Output Argument: 1266 . q - A PetscQuadrature object 1267 1268 Level: intermediate 1269 1270 .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval() 1271 @*/ 1272 PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 1273 { 1274 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c; 1275 PetscReal *x, *w, *xw, *ww; 1276 PetscErrorCode ierr; 1277 1278 PetscFunctionBegin; 1279 ierr = PetscMalloc1(totpoints*dim,&x);CHKERRQ(ierr); 1280 ierr = PetscMalloc1(totpoints*Nc,&w);CHKERRQ(ierr); 1281 /* Set up the Golub-Welsch system */ 1282 switch (dim) { 1283 case 0: 1284 ierr = PetscFree(x);CHKERRQ(ierr); 1285 ierr = PetscFree(w);CHKERRQ(ierr); 1286 ierr = PetscMalloc1(1, &x);CHKERRQ(ierr); 1287 ierr = PetscMalloc1(Nc, &w);CHKERRQ(ierr); 1288 x[0] = 0.0; 1289 for (c = 0; c < Nc; ++c) w[c] = 1.0; 1290 break; 1291 case 1: 1292 ierr = PetscMalloc1(npoints,&ww);CHKERRQ(ierr); 1293 ierr = PetscDTGaussQuadrature(npoints, a, b, x, ww);CHKERRQ(ierr); 1294 for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i]; 1295 ierr = PetscFree(ww);CHKERRQ(ierr); 1296 break; 1297 case 2: 1298 ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr); 1299 ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr); 1300 for (i = 0; i < npoints; ++i) { 1301 for (j = 0; j < npoints; ++j) { 1302 x[(i*npoints+j)*dim+0] = xw[i]; 1303 x[(i*npoints+j)*dim+1] = xw[j]; 1304 for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j]; 1305 } 1306 } 1307 ierr = PetscFree2(xw,ww);CHKERRQ(ierr); 1308 break; 1309 case 3: 1310 ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr); 1311 ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr); 1312 for (i = 0; i < npoints; ++i) { 1313 for (j = 0; j < npoints; ++j) { 1314 for (k = 0; k < npoints; ++k) { 1315 x[((i*npoints+j)*npoints+k)*dim+0] = xw[i]; 1316 x[((i*npoints+j)*npoints+k)*dim+1] = xw[j]; 1317 x[((i*npoints+j)*npoints+k)*dim+2] = xw[k]; 1318 for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k]; 1319 } 1320 } 1321 } 1322 ierr = PetscFree2(xw,ww);CHKERRQ(ierr); 1323 break; 1324 default: 1325 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim); 1326 } 1327 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 1328 ierr = PetscQuadratureSetOrder(*q, 2*npoints-1);CHKERRQ(ierr); 1329 ierr = PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);CHKERRQ(ierr); 1330 ierr = PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");CHKERRQ(ierr); 1331 PetscFunctionReturn(0); 1332 } 1333 1334 /* Maps from [-1,1]^2 to the (-1,1) reference triangle */ 1335 PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta) 1336 { 1337 PetscFunctionBegin; 1338 *xi = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0; 1339 *eta = 0.5 * (1.0 + y) * (1.0 - z) - 1.0; 1340 *zeta = z; 1341 PetscFunctionReturn(0); 1342 } 1343 1344 1345 /*@ 1346 PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex 1347 1348 Not Collective 1349 1350 Input Arguments: 1351 + dim - The simplex dimension 1352 . Nc - The number of components 1353 . npoints - The number of points in one dimension 1354 . a - left end of interval (often-1) 1355 - b - right end of interval (often +1) 1356 1357 Output Argument: 1358 . q - A PetscQuadrature object 1359 1360 Level: intermediate 1361 1362 References: 1363 . 1. - Karniadakis and Sherwin. FIAT 1364 1365 Note: For dim == 1, this is Gauss-Legendre quadrature 1366 1367 .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature() 1368 @*/ 1369 PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 1370 { 1371 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints; 1372 PetscReal *px, *wx, *py, *wy, *pz, *wz, *x, *w; 1373 PetscInt i, j, k, c; PetscErrorCode ierr; 1374 1375 PetscFunctionBegin; 1376 if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now"); 1377 ierr = PetscMalloc1(totpoints*dim, &x);CHKERRQ(ierr); 1378 ierr = PetscMalloc1(totpoints*Nc, &w);CHKERRQ(ierr); 1379 switch (dim) { 1380 case 0: 1381 ierr = PetscFree(x);CHKERRQ(ierr); 1382 ierr = PetscFree(w);CHKERRQ(ierr); 1383 ierr = PetscMalloc1(1, &x);CHKERRQ(ierr); 1384 ierr = PetscMalloc1(Nc, &w);CHKERRQ(ierr); 1385 x[0] = 0.0; 1386 for (c = 0; c < Nc; ++c) w[c] = 1.0; 1387 break; 1388 case 1: 1389 ierr = PetscMalloc1(npoints,&wx);CHKERRQ(ierr); 1390 ierr = PetscDTGaussJacobiQuadrature(npoints, -1., 1., 0.0, 0.0, x, wx);CHKERRQ(ierr); 1391 for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i]; 1392 ierr = PetscFree(wx);CHKERRQ(ierr); 1393 break; 1394 case 2: 1395 ierr = PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);CHKERRQ(ierr); 1396 ierr = PetscDTGaussJacobiQuadrature(npoints, -1., 1., 0.0, 0.0, px, wx);CHKERRQ(ierr); 1397 ierr = PetscDTGaussJacobiQuadrature(npoints, -1., 1., 1.0, 0.0, py, wy);CHKERRQ(ierr); 1398 for (i = 0; i < npoints; ++i) { 1399 for (j = 0; j < npoints; ++j) { 1400 ierr = PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);CHKERRQ(ierr); 1401 for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j]; 1402 } 1403 } 1404 ierr = PetscFree4(px,wx,py,wy);CHKERRQ(ierr); 1405 break; 1406 case 3: 1407 ierr = PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);CHKERRQ(ierr); 1408 ierr = PetscDTGaussJacobiQuadrature(npoints, -1., 1., 0.0, 0.0, px, wx);CHKERRQ(ierr); 1409 ierr = PetscDTGaussJacobiQuadrature(npoints, -1., 1., 1.0, 0.0, py, wy);CHKERRQ(ierr); 1410 ierr = PetscDTGaussJacobiQuadrature(npoints, -1., 1., 2.0, 0.0, pz, wz);CHKERRQ(ierr); 1411 for (i = 0; i < npoints; ++i) { 1412 for (j = 0; j < npoints; ++j) { 1413 for (k = 0; k < npoints; ++k) { 1414 ierr = PetscDTMapCubeToTetrahedron_Internal(px[i], py[j], pz[k], &x[((i*npoints+j)*npoints+k)*3+0], &x[((i*npoints+j)*npoints+k)*3+1], &x[((i*npoints+j)*npoints+k)*3+2]);CHKERRQ(ierr); 1415 for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k]; 1416 } 1417 } 1418 } 1419 ierr = PetscFree6(px,wx,py,wy,pz,wz);CHKERRQ(ierr); 1420 break; 1421 default: 1422 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim); 1423 } 1424 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 1425 ierr = PetscQuadratureSetOrder(*q, 2*npoints-1);CHKERRQ(ierr); 1426 ierr = PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);CHKERRQ(ierr); 1427 ierr = PetscObjectChangeTypeName((PetscObject)*q,"GaussJacobi");CHKERRQ(ierr); 1428 PetscFunctionReturn(0); 1429 } 1430 1431 /*@ 1432 PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell 1433 1434 Not Collective 1435 1436 Input Arguments: 1437 + dim - The cell dimension 1438 . level - The number of points in one dimension, 2^l 1439 . a - left end of interval (often-1) 1440 - b - right end of interval (often +1) 1441 1442 Output Argument: 1443 . q - A PetscQuadrature object 1444 1445 Level: intermediate 1446 1447 .seealso: PetscDTGaussTensorQuadrature() 1448 @*/ 1449 PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q) 1450 { 1451 const PetscInt p = 16; /* Digits of precision in the evaluation */ 1452 const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */ 1453 const PetscReal beta = (b+a)/2.; /* Center of the integration interval */ 1454 const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */ 1455 PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */ 1456 PetscReal wk = 0.5*PETSC_PI; /* Quadrature weight at x_k */ 1457 PetscReal *x, *w; 1458 PetscInt K, k, npoints; 1459 PetscErrorCode ierr; 1460 1461 PetscFunctionBegin; 1462 if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim); 1463 if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits"); 1464 /* Find K such that the weights are < 32 digits of precision */ 1465 for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) { 1466 wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h))); 1467 } 1468 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 1469 ierr = PetscQuadratureSetOrder(*q, 2*K+1);CHKERRQ(ierr); 1470 npoints = 2*K-1; 1471 ierr = PetscMalloc1(npoints*dim, &x);CHKERRQ(ierr); 1472 ierr = PetscMalloc1(npoints, &w);CHKERRQ(ierr); 1473 /* Center term */ 1474 x[0] = beta; 1475 w[0] = 0.5*alpha*PETSC_PI; 1476 for (k = 1; k < K; ++k) { 1477 wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1478 xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h)); 1479 x[2*k-1] = -alpha*xk+beta; 1480 w[2*k-1] = wk; 1481 x[2*k+0] = alpha*xk+beta; 1482 w[2*k+0] = wk; 1483 } 1484 ierr = PetscQuadratureSetData(*q, dim, 1, npoints, x, w);CHKERRQ(ierr); 1485 PetscFunctionReturn(0); 1486 } 1487 1488 PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 1489 { 1490 const PetscInt p = 16; /* Digits of precision in the evaluation */ 1491 const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */ 1492 const PetscReal beta = (b+a)/2.; /* Center of the integration interval */ 1493 PetscReal h = 1.0; /* Step size, length between x_k */ 1494 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 1495 PetscReal osum = 0.0; /* Integral on last level */ 1496 PetscReal psum = 0.0; /* Integral on the level before the last level */ 1497 PetscReal sum; /* Integral on current level */ 1498 PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 1499 PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 1500 PetscReal wk; /* Quadrature weight at x_k */ 1501 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 1502 PetscInt d; /* Digits of precision in the integral */ 1503 1504 PetscFunctionBegin; 1505 if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 1506 /* Center term */ 1507 func(beta, &lval); 1508 sum = 0.5*alpha*PETSC_PI*lval; 1509 /* */ 1510 do { 1511 PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4; 1512 PetscInt k = 1; 1513 1514 ++l; 1515 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */ 1516 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 1517 psum = osum; 1518 osum = sum; 1519 h *= 0.5; 1520 sum *= 0.5; 1521 do { 1522 wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1523 yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1524 lx = -alpha*(1.0 - yk)+beta; 1525 rx = alpha*(1.0 - yk)+beta; 1526 func(lx, &lval); 1527 func(rx, &rval); 1528 lterm = alpha*wk*lval; 1529 maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm); 1530 sum += lterm; 1531 rterm = alpha*wk*rval; 1532 maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm); 1533 sum += rterm; 1534 ++k; 1535 /* Only need to evaluate every other point on refined levels */ 1536 if (l != 1) ++k; 1537 } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */ 1538 1539 d1 = PetscLog10Real(PetscAbsReal(sum - osum)); 1540 d2 = PetscLog10Real(PetscAbsReal(sum - psum)); 1541 d3 = PetscLog10Real(maxTerm) - p; 1542 if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0; 1543 else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm))); 1544 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4))); 1545 } while (d < digits && l < 12); 1546 *sol = sum; 1547 1548 PetscFunctionReturn(0); 1549 } 1550 1551 #if defined(PETSC_HAVE_MPFR) 1552 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 1553 { 1554 const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */ 1555 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 1556 mpfr_t alpha; /* Half-width of the integration interval */ 1557 mpfr_t beta; /* Center of the integration interval */ 1558 mpfr_t h; /* Step size, length between x_k */ 1559 mpfr_t osum; /* Integral on last level */ 1560 mpfr_t psum; /* Integral on the level before the last level */ 1561 mpfr_t sum; /* Integral on current level */ 1562 mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 1563 mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 1564 mpfr_t wk; /* Quadrature weight at x_k */ 1565 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 1566 PetscInt d; /* Digits of precision in the integral */ 1567 mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp; 1568 1569 PetscFunctionBegin; 1570 if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 1571 /* Create high precision storage */ 1572 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); 1573 /* Initialization */ 1574 mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN); 1575 mpfr_set_d(beta, 0.5*(b+a), MPFR_RNDN); 1576 mpfr_set_d(osum, 0.0, MPFR_RNDN); 1577 mpfr_set_d(psum, 0.0, MPFR_RNDN); 1578 mpfr_set_d(h, 1.0, MPFR_RNDN); 1579 mpfr_const_pi(pi2, MPFR_RNDN); 1580 mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN); 1581 /* Center term */ 1582 func(0.5*(b+a), &lval); 1583 mpfr_set(sum, pi2, MPFR_RNDN); 1584 mpfr_mul(sum, sum, alpha, MPFR_RNDN); 1585 mpfr_mul_d(sum, sum, lval, MPFR_RNDN); 1586 /* */ 1587 do { 1588 PetscReal d1, d2, d3, d4; 1589 PetscInt k = 1; 1590 1591 ++l; 1592 mpfr_set_d(maxTerm, 0.0, MPFR_RNDN); 1593 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */ 1594 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 1595 mpfr_set(psum, osum, MPFR_RNDN); 1596 mpfr_set(osum, sum, MPFR_RNDN); 1597 mpfr_mul_d(h, h, 0.5, MPFR_RNDN); 1598 mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN); 1599 do { 1600 mpfr_set_si(kh, k, MPFR_RNDN); 1601 mpfr_mul(kh, kh, h, MPFR_RNDN); 1602 /* Weight */ 1603 mpfr_set(wk, h, MPFR_RNDN); 1604 mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN); 1605 mpfr_mul(msinh, msinh, pi2, MPFR_RNDN); 1606 mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN); 1607 mpfr_cosh(tmp, msinh, MPFR_RNDN); 1608 mpfr_sqr(tmp, tmp, MPFR_RNDN); 1609 mpfr_mul(wk, wk, mcosh, MPFR_RNDN); 1610 mpfr_div(wk, wk, tmp, MPFR_RNDN); 1611 /* Abscissa */ 1612 mpfr_set_d(yk, 1.0, MPFR_RNDZ); 1613 mpfr_cosh(tmp, msinh, MPFR_RNDN); 1614 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 1615 mpfr_exp(tmp, msinh, MPFR_RNDN); 1616 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 1617 /* Quadrature points */ 1618 mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ); 1619 mpfr_mul(lx, lx, alpha, MPFR_RNDU); 1620 mpfr_add(lx, lx, beta, MPFR_RNDU); 1621 mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ); 1622 mpfr_mul(rx, rx, alpha, MPFR_RNDD); 1623 mpfr_add(rx, rx, beta, MPFR_RNDD); 1624 /* Evaluation */ 1625 func(mpfr_get_d(lx, MPFR_RNDU), &lval); 1626 func(mpfr_get_d(rx, MPFR_RNDD), &rval); 1627 /* Update */ 1628 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 1629 mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN); 1630 mpfr_add(sum, sum, tmp, MPFR_RNDN); 1631 mpfr_abs(tmp, tmp, MPFR_RNDN); 1632 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 1633 mpfr_set(curTerm, tmp, MPFR_RNDN); 1634 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 1635 mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN); 1636 mpfr_add(sum, sum, tmp, MPFR_RNDN); 1637 mpfr_abs(tmp, tmp, MPFR_RNDN); 1638 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 1639 mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN); 1640 ++k; 1641 /* Only need to evaluate every other point on refined levels */ 1642 if (l != 1) ++k; 1643 mpfr_log10(tmp, wk, MPFR_RNDN); 1644 mpfr_abs(tmp, tmp, MPFR_RNDN); 1645 } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */ 1646 mpfr_sub(tmp, sum, osum, MPFR_RNDN); 1647 mpfr_abs(tmp, tmp, MPFR_RNDN); 1648 mpfr_log10(tmp, tmp, MPFR_RNDN); 1649 d1 = mpfr_get_d(tmp, MPFR_RNDN); 1650 mpfr_sub(tmp, sum, psum, MPFR_RNDN); 1651 mpfr_abs(tmp, tmp, MPFR_RNDN); 1652 mpfr_log10(tmp, tmp, MPFR_RNDN); 1653 d2 = mpfr_get_d(tmp, MPFR_RNDN); 1654 mpfr_log10(tmp, maxTerm, MPFR_RNDN); 1655 d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits; 1656 mpfr_log10(tmp, curTerm, MPFR_RNDN); 1657 d4 = mpfr_get_d(tmp, MPFR_RNDN); 1658 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4))); 1659 } while (d < digits && l < 8); 1660 *sol = mpfr_get_d(sum, MPFR_RNDN); 1661 /* Cleanup */ 1662 mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL); 1663 PetscFunctionReturn(0); 1664 } 1665 #else 1666 1667 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 1668 { 1669 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp"); 1670 } 1671 #endif 1672 1673 /* Overwrites A. Can only handle full-rank problems with m>=n 1674 * A in column-major format 1675 * Ainv in row-major format 1676 * tau has length m 1677 * worksize must be >= max(1,n) 1678 */ 1679 static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work) 1680 { 1681 PetscErrorCode ierr; 1682 PetscBLASInt M,N,K,lda,ldb,ldwork,info; 1683 PetscScalar *A,*Ainv,*R,*Q,Alpha; 1684 1685 PetscFunctionBegin; 1686 #if defined(PETSC_USE_COMPLEX) 1687 { 1688 PetscInt i,j; 1689 ierr = PetscMalloc2(m*n,&A,m*n,&Ainv);CHKERRQ(ierr); 1690 for (j=0; j<n; j++) { 1691 for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j]; 1692 } 1693 mstride = m; 1694 } 1695 #else 1696 A = A_in; 1697 Ainv = Ainv_out; 1698 #endif 1699 1700 ierr = PetscBLASIntCast(m,&M);CHKERRQ(ierr); 1701 ierr = PetscBLASIntCast(n,&N);CHKERRQ(ierr); 1702 ierr = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr); 1703 ierr = PetscBLASIntCast(worksize,&ldwork);CHKERRQ(ierr); 1704 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1705 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info)); 1706 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1707 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error"); 1708 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */ 1709 1710 /* Extract an explicit representation of Q */ 1711 Q = Ainv; 1712 ierr = PetscArraycpy(Q,A,mstride*n);CHKERRQ(ierr); 1713 K = N; /* full rank */ 1714 PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info)); 1715 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error"); 1716 1717 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */ 1718 Alpha = 1.0; 1719 ldb = lda; 1720 PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb)); 1721 /* Ainv is Q, overwritten with inverse */ 1722 1723 #if defined(PETSC_USE_COMPLEX) 1724 { 1725 PetscInt i; 1726 for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]); 1727 ierr = PetscFree2(A,Ainv);CHKERRQ(ierr); 1728 } 1729 #endif 1730 PetscFunctionReturn(0); 1731 } 1732 1733 /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */ 1734 static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B) 1735 { 1736 PetscErrorCode ierr; 1737 PetscReal *Bv; 1738 PetscInt i,j; 1739 1740 PetscFunctionBegin; 1741 ierr = PetscMalloc1((ninterval+1)*ndegree,&Bv);CHKERRQ(ierr); 1742 /* Point evaluation of L_p on all the source vertices */ 1743 ierr = PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);CHKERRQ(ierr); 1744 /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */ 1745 for (i=0; i<ninterval; i++) { 1746 for (j=0; j<ndegree; j++) { 1747 if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j]; 1748 else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j]; 1749 } 1750 } 1751 ierr = PetscFree(Bv);CHKERRQ(ierr); 1752 PetscFunctionReturn(0); 1753 } 1754 1755 /*@ 1756 PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals 1757 1758 Not Collective 1759 1760 Input Arguments: 1761 + degree - degree of reconstruction polynomial 1762 . nsource - number of source intervals 1763 . sourcex - sorted coordinates of source cell boundaries (length nsource+1) 1764 . ntarget - number of target intervals 1765 - targetx - sorted coordinates of target cell boundaries (length ntarget+1) 1766 1767 Output Arguments: 1768 . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s] 1769 1770 Level: advanced 1771 1772 .seealso: PetscDTLegendreEval() 1773 @*/ 1774 PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R) 1775 { 1776 PetscErrorCode ierr; 1777 PetscInt i,j,k,*bdegrees,worksize; 1778 PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget; 1779 PetscScalar *tau,*work; 1780 1781 PetscFunctionBegin; 1782 PetscValidRealPointer(sourcex,3); 1783 PetscValidRealPointer(targetx,5); 1784 PetscValidRealPointer(R,6); 1785 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); 1786 #if defined(PETSC_USE_DEBUG) 1787 for (i=0; i<nsource; i++) { 1788 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]); 1789 } 1790 for (i=0; i<ntarget; i++) { 1791 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]); 1792 } 1793 #endif 1794 xmin = PetscMin(sourcex[0],targetx[0]); 1795 xmax = PetscMax(sourcex[nsource],targetx[ntarget]); 1796 center = (xmin + xmax)/2; 1797 hscale = (xmax - xmin)/2; 1798 worksize = nsource; 1799 ierr = PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);CHKERRQ(ierr); 1800 ierr = PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);CHKERRQ(ierr); 1801 for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale; 1802 for (i=0; i<=degree; i++) bdegrees[i] = i+1; 1803 ierr = PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);CHKERRQ(ierr); 1804 ierr = PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);CHKERRQ(ierr); 1805 for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale; 1806 ierr = PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);CHKERRQ(ierr); 1807 for (i=0; i<ntarget; i++) { 1808 PetscReal rowsum = 0; 1809 for (j=0; j<nsource; j++) { 1810 PetscReal sum = 0; 1811 for (k=0; k<degree+1; k++) { 1812 sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j]; 1813 } 1814 R[i*nsource+j] = sum; 1815 rowsum += sum; 1816 } 1817 for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */ 1818 } 1819 ierr = PetscFree4(bdegrees,sourcey,Bsource,work);CHKERRQ(ierr); 1820 ierr = PetscFree4(tau,Bsinv,targety,Btarget);CHKERRQ(ierr); 1821 PetscFunctionReturn(0); 1822 } 1823 1824 /*@C 1825 PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points 1826 1827 Not Collective 1828 1829 Input Parameter: 1830 + n - the number of GLL nodes 1831 . nodes - the GLL nodes 1832 . weights - the GLL weights 1833 - f - the function values at the nodes 1834 1835 Output Parameter: 1836 . in - the value of the integral 1837 1838 Level: beginner 1839 1840 .seealso: PetscDTGaussLobattoLegendreQuadrature() 1841 1842 @*/ 1843 PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal *nodes,PetscReal *weights,const PetscReal *f,PetscReal *in) 1844 { 1845 PetscInt i; 1846 1847 PetscFunctionBegin; 1848 *in = 0.; 1849 for (i=0; i<n; i++) { 1850 *in += f[i]*f[i]*weights[i]; 1851 } 1852 PetscFunctionReturn(0); 1853 } 1854 1855 /*@C 1856 PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element 1857 1858 Not Collective 1859 1860 Input Parameter: 1861 + n - the number of GLL nodes 1862 . nodes - the GLL nodes 1863 - weights - the GLL weights 1864 1865 Output Parameter: 1866 . A - the stiffness element 1867 1868 Level: beginner 1869 1870 Notes: 1871 Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy() 1872 1873 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) 1874 1875 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy() 1876 1877 @*/ 1878 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1879 { 1880 PetscReal **A; 1881 PetscErrorCode ierr; 1882 const PetscReal *gllnodes = nodes; 1883 const PetscInt p = n-1; 1884 PetscReal z0,z1,z2 = -1,x,Lpj,Lpr; 1885 PetscInt i,j,nn,r; 1886 1887 PetscFunctionBegin; 1888 ierr = PetscMalloc1(n,&A);CHKERRQ(ierr); 1889 ierr = PetscMalloc1(n*n,&A[0]);CHKERRQ(ierr); 1890 for (i=1; i<n; i++) A[i] = A[i-1]+n; 1891 1892 for (j=1; j<p; j++) { 1893 x = gllnodes[j]; 1894 z0 = 1.; 1895 z1 = x; 1896 for (nn=1; nn<p; nn++) { 1897 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1898 z0 = z1; 1899 z1 = z2; 1900 } 1901 Lpj=z2; 1902 for (r=1; r<p; r++) { 1903 if (r == j) { 1904 A[j][j]=2./(3.*(1.-gllnodes[j]*gllnodes[j])*Lpj*Lpj); 1905 } else { 1906 x = gllnodes[r]; 1907 z0 = 1.; 1908 z1 = x; 1909 for (nn=1; nn<p; nn++) { 1910 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1911 z0 = z1; 1912 z1 = z2; 1913 } 1914 Lpr = z2; 1915 A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(gllnodes[j]-gllnodes[r])*(gllnodes[j]-gllnodes[r])); 1916 } 1917 } 1918 } 1919 for (j=1; j<p+1; j++) { 1920 x = gllnodes[j]; 1921 z0 = 1.; 1922 z1 = x; 1923 for (nn=1; nn<p; nn++) { 1924 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1925 z0 = z1; 1926 z1 = z2; 1927 } 1928 Lpj = z2; 1929 A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+gllnodes[j])*(1.+gllnodes[j])); 1930 A[0][j] = A[j][0]; 1931 } 1932 for (j=0; j<p; j++) { 1933 x = gllnodes[j]; 1934 z0 = 1.; 1935 z1 = x; 1936 for (nn=1; nn<p; nn++) { 1937 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1938 z0 = z1; 1939 z1 = z2; 1940 } 1941 Lpj=z2; 1942 1943 A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-gllnodes[j])*(1.-gllnodes[j])); 1944 A[j][p] = A[p][j]; 1945 } 1946 A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.; 1947 A[p][p]=A[0][0]; 1948 *AA = A; 1949 PetscFunctionReturn(0); 1950 } 1951 1952 /*@C 1953 PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element 1954 1955 Not Collective 1956 1957 Input Parameter: 1958 + n - the number of GLL nodes 1959 . nodes - the GLL nodes 1960 . weights - the GLL weightss 1961 - A - the stiffness element 1962 1963 Level: beginner 1964 1965 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate() 1966 1967 @*/ 1968 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1969 { 1970 PetscErrorCode ierr; 1971 1972 PetscFunctionBegin; 1973 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 1974 ierr = PetscFree(*AA);CHKERRQ(ierr); 1975 *AA = NULL; 1976 PetscFunctionReturn(0); 1977 } 1978 1979 /*@C 1980 PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element 1981 1982 Not Collective 1983 1984 Input Parameter: 1985 + n - the number of GLL nodes 1986 . nodes - the GLL nodes 1987 . weights - the GLL weights 1988 1989 Output Parameter: 1990 . AA - the stiffness element 1991 - AAT - the transpose of AA (pass in NULL if you do not need this array) 1992 1993 Level: beginner 1994 1995 Notes: 1996 Destroy this with PetscGaussLobattoLegendreElementGradientDestroy() 1997 1998 You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented 1999 2000 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy() 2001 2002 @*/ 2003 PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT) 2004 { 2005 PetscReal **A, **AT = NULL; 2006 PetscErrorCode ierr; 2007 const PetscReal *gllnodes = nodes; 2008 const PetscInt p = n-1; 2009 PetscReal Li, Lj,d0; 2010 PetscInt i,j; 2011 2012 PetscFunctionBegin; 2013 ierr = PetscMalloc1(n,&A);CHKERRQ(ierr); 2014 ierr = PetscMalloc1(n*n,&A[0]);CHKERRQ(ierr); 2015 for (i=1; i<n; i++) A[i] = A[i-1]+n; 2016 2017 if (AAT) { 2018 ierr = PetscMalloc1(n,&AT);CHKERRQ(ierr); 2019 ierr = PetscMalloc1(n*n,&AT[0]);CHKERRQ(ierr); 2020 for (i=1; i<n; i++) AT[i] = AT[i-1]+n; 2021 } 2022 2023 if (n==1) {A[0][0] = 0.;} 2024 d0 = (PetscReal)p*((PetscReal)p+1.)/4.; 2025 for (i=0; i<n; i++) { 2026 for (j=0; j<n; j++) { 2027 A[i][j] = 0.; 2028 ierr = PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li);CHKERRQ(ierr); 2029 ierr = PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj);CHKERRQ(ierr); 2030 if (i!=j) A[i][j] = Li/(Lj*(gllnodes[i]-gllnodes[j])); 2031 if ((j==i) && (i==0)) A[i][j] = -d0; 2032 if (j==i && i==p) A[i][j] = d0; 2033 if (AT) AT[j][i] = A[i][j]; 2034 } 2035 } 2036 if (AAT) *AAT = AT; 2037 *AA = A; 2038 PetscFunctionReturn(0); 2039 } 2040 2041 /*@C 2042 PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGaussLobattoLegendreElementGradientCreate() 2043 2044 Not Collective 2045 2046 Input Parameter: 2047 + n - the number of GLL nodes 2048 . nodes - the GLL nodes 2049 . weights - the GLL weights 2050 . AA - the stiffness element 2051 - AAT - the transpose of the element 2052 2053 Level: beginner 2054 2055 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionCreate() 2056 2057 @*/ 2058 PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT) 2059 { 2060 PetscErrorCode ierr; 2061 2062 PetscFunctionBegin; 2063 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 2064 ierr = PetscFree(*AA);CHKERRQ(ierr); 2065 *AA = NULL; 2066 if (*AAT) { 2067 ierr = PetscFree((*AAT)[0]);CHKERRQ(ierr); 2068 ierr = PetscFree(*AAT);CHKERRQ(ierr); 2069 *AAT = NULL; 2070 } 2071 PetscFunctionReturn(0); 2072 } 2073 2074 /*@C 2075 PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element 2076 2077 Not Collective 2078 2079 Input Parameter: 2080 + n - the number of GLL nodes 2081 . nodes - the GLL nodes 2082 - weights - the GLL weightss 2083 2084 Output Parameter: 2085 . AA - the stiffness element 2086 2087 Level: beginner 2088 2089 Notes: 2090 Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy() 2091 2092 This is the same as the Gradient operator multiplied by the diagonal mass matrix 2093 2094 You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented 2095 2096 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionDestroy() 2097 2098 @*/ 2099 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 2100 { 2101 PetscReal **D; 2102 PetscErrorCode ierr; 2103 const PetscReal *gllweights = weights; 2104 const PetscInt glln = n; 2105 PetscInt i,j; 2106 2107 PetscFunctionBegin; 2108 ierr = PetscGaussLobattoLegendreElementGradientCreate(n,nodes,weights,&D,NULL);CHKERRQ(ierr); 2109 for (i=0; i<glln; i++){ 2110 for (j=0; j<glln; j++) { 2111 D[i][j] = gllweights[i]*D[i][j]; 2112 } 2113 } 2114 *AA = D; 2115 PetscFunctionReturn(0); 2116 } 2117 2118 /*@C 2119 PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element 2120 2121 Not Collective 2122 2123 Input Parameter: 2124 + n - the number of GLL nodes 2125 . nodes - the GLL nodes 2126 . weights - the GLL weights 2127 - A - advection 2128 2129 Level: beginner 2130 2131 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementAdvectionCreate() 2132 2133 @*/ 2134 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 2135 { 2136 PetscErrorCode ierr; 2137 2138 PetscFunctionBegin; 2139 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 2140 ierr = PetscFree(*AA);CHKERRQ(ierr); 2141 *AA = NULL; 2142 PetscFunctionReturn(0); 2143 } 2144 2145 PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 2146 { 2147 PetscReal **A; 2148 PetscErrorCode ierr; 2149 const PetscReal *gllweights = weights; 2150 const PetscInt glln = n; 2151 PetscInt i,j; 2152 2153 PetscFunctionBegin; 2154 ierr = PetscMalloc1(glln,&A);CHKERRQ(ierr); 2155 ierr = PetscMalloc1(glln*glln,&A[0]);CHKERRQ(ierr); 2156 for (i=1; i<glln; i++) A[i] = A[i-1]+glln; 2157 if (glln==1) {A[0][0] = 0.;} 2158 for (i=0; i<glln; i++) { 2159 for (j=0; j<glln; j++) { 2160 A[i][j] = 0.; 2161 if (j==i) A[i][j] = gllweights[i]; 2162 } 2163 } 2164 *AA = A; 2165 PetscFunctionReturn(0); 2166 } 2167 2168 PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 2169 { 2170 PetscErrorCode ierr; 2171 2172 PetscFunctionBegin; 2173 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 2174 ierr = PetscFree(*AA);CHKERRQ(ierr); 2175 *AA = NULL; 2176 PetscFunctionReturn(0); 2177 } 2178 2179 /*@ 2180 PetscDTIndexToBary - convert an index into a barycentric coordinate. 2181 2182 Input Parameters: 2183 + 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) 2184 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 2185 - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum) 2186 2187 Output Parameter: 2188 . coord - will be filled with the barycentric coordinate 2189 2190 Level: beginner 2191 2192 Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the 2193 least significant and the last index is the most significant. 2194 2195 .seealso: PetscDTBaryToIndex 2196 @*/ 2197 PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[]) 2198 { 2199 PetscInt c, d, s, total, subtotal, nexttotal; 2200 2201 PetscFunctionBeginHot; 2202 if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 2203 if (index < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative"); 2204 if (!len) { 2205 if (!sum && !index) PetscFunctionReturn(0); 2206 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 2207 } 2208 for (c = 1, total = 1; c <= len; c++) { 2209 /* total is the number of ways to have a tuple of length c with sum */ 2210 if (index < total) break; 2211 total = (total * (sum + c)) / c; 2212 } 2213 if (c > len) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range"); 2214 for (d = c; d < len; d++) coord[d] = 0; 2215 for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) { 2216 /* subtotal is the number of ways to have a tuple of length c with sum s */ 2217 /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */ 2218 if ((index + subtotal) >= total) { 2219 coord[--c] = sum - s; 2220 index -= (total - subtotal); 2221 sum = s; 2222 total = nexttotal; 2223 subtotal = 1; 2224 nexttotal = 1; 2225 s = 0; 2226 } else { 2227 subtotal = (subtotal * (c + s)) / (s + 1); 2228 nexttotal = (nexttotal * (c - 1 + s)) / (s + 1); 2229 s++; 2230 } 2231 } 2232 PetscFunctionReturn(0); 2233 } 2234 2235 /*@ 2236 PetscDTBaryToIndex - convert a barycentric coordinate to an index 2237 2238 Input Parameters: 2239 + 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) 2240 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to 2241 - coord - a barycentric coordinate with the given length and sum 2242 2243 Output Parameter: 2244 . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum) 2245 2246 Level: beginner 2247 2248 Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the 2249 least significant and the last index is the most significant. 2250 2251 .seealso: PetscDTIndexToBary 2252 @*/ 2253 PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index) 2254 { 2255 PetscInt c; 2256 PetscInt i; 2257 PetscInt total; 2258 2259 PetscFunctionBeginHot; 2260 if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative"); 2261 if (!len) { 2262 if (!sum) { 2263 *index = 0; 2264 PetscFunctionReturn(0); 2265 } 2266 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate"); 2267 } 2268 for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c; 2269 i = total - 1; 2270 c = len - 1; 2271 sum -= coord[c]; 2272 while (sum > 0) { 2273 PetscInt subtotal; 2274 PetscInt s; 2275 2276 for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s; 2277 i -= subtotal; 2278 sum -= coord[--c]; 2279 } 2280 *index = i; 2281 PetscFunctionReturn(0); 2282 } 2283