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