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 static PetscBool GaussCite = PETSC_FALSE; 16 const char GaussCitation[] = "@article{GolubWelsch1969,\n" 17 " author = {Golub and Welsch},\n" 18 " title = {Calculation of Quadrature Rules},\n" 19 " journal = {Math. Comp.},\n" 20 " volume = {23},\n" 21 " number = {106},\n" 22 " pages = {221--230},\n" 23 " year = {1969}\n}\n"; 24 25 26 PetscClassId PETSCQUADRATURE_CLASSID = 0; 27 28 /*@ 29 PetscQuadratureCreate - Create a PetscQuadrature object 30 31 Collective 32 33 Input Parameter: 34 . comm - The communicator for the PetscQuadrature object 35 36 Output Parameter: 37 . q - The PetscQuadrature object 38 39 Level: beginner 40 41 .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData() 42 @*/ 43 PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q) 44 { 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 PetscValidPointer(q, 2); 49 ierr = DMInitializePackage();CHKERRQ(ierr); 50 ierr = PetscHeaderCreate(*q,PETSCQUADRATURE_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);CHKERRQ(ierr); 51 (*q)->dim = -1; 52 (*q)->Nc = 1; 53 (*q)->order = -1; 54 (*q)->numPoints = 0; 55 (*q)->points = NULL; 56 (*q)->weights = NULL; 57 PetscFunctionReturn(0); 58 } 59 60 /*@ 61 PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object 62 63 Collective on q 64 65 Input Parameter: 66 . q - The PetscQuadrature object 67 68 Output Parameter: 69 . r - The new PetscQuadrature object 70 71 Level: beginner 72 73 .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData() 74 @*/ 75 PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r) 76 { 77 PetscInt order, dim, Nc, Nq; 78 const PetscReal *points, *weights; 79 PetscReal *p, *w; 80 PetscErrorCode ierr; 81 82 PetscFunctionBegin; 83 PetscValidPointer(q, 2); 84 ierr = PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);CHKERRQ(ierr); 85 ierr = PetscQuadratureGetOrder(q, &order);CHKERRQ(ierr); 86 ierr = PetscQuadratureSetOrder(*r, order);CHKERRQ(ierr); 87 ierr = PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);CHKERRQ(ierr); 88 ierr = PetscMalloc1(Nq*dim, &p);CHKERRQ(ierr); 89 ierr = PetscMalloc1(Nq*Nc, &w);CHKERRQ(ierr); 90 ierr = PetscArraycpy(p, points, Nq*dim);CHKERRQ(ierr); 91 ierr = PetscArraycpy(w, weights, Nc * Nq);CHKERRQ(ierr); 92 ierr = PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 /*@ 97 PetscQuadratureDestroy - Destroys a PetscQuadrature object 98 99 Collective on q 100 101 Input Parameter: 102 . q - The PetscQuadrature object 103 104 Level: beginner 105 106 .seealso: PetscQuadratureCreate(), PetscQuadratureGetData() 107 @*/ 108 PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q) 109 { 110 PetscErrorCode ierr; 111 112 PetscFunctionBegin; 113 if (!*q) PetscFunctionReturn(0); 114 PetscValidHeaderSpecific((*q),PETSCQUADRATURE_CLASSID,1); 115 if (--((PetscObject)(*q))->refct > 0) { 116 *q = NULL; 117 PetscFunctionReturn(0); 118 } 119 ierr = PetscFree((*q)->points);CHKERRQ(ierr); 120 ierr = PetscFree((*q)->weights);CHKERRQ(ierr); 121 ierr = PetscHeaderDestroy(q);CHKERRQ(ierr); 122 PetscFunctionReturn(0); 123 } 124 125 /*@ 126 PetscQuadratureGetOrder - Return the order of the method 127 128 Not collective 129 130 Input Parameter: 131 . q - The PetscQuadrature object 132 133 Output Parameter: 134 . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated 135 136 Level: intermediate 137 138 .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData() 139 @*/ 140 PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order) 141 { 142 PetscFunctionBegin; 143 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 144 PetscValidPointer(order, 2); 145 *order = q->order; 146 PetscFunctionReturn(0); 147 } 148 149 /*@ 150 PetscQuadratureSetOrder - Return the order of the method 151 152 Not collective 153 154 Input Parameters: 155 + q - The PetscQuadrature object 156 - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated 157 158 Level: intermediate 159 160 .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData() 161 @*/ 162 PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order) 163 { 164 PetscFunctionBegin; 165 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 166 q->order = order; 167 PetscFunctionReturn(0); 168 } 169 170 /*@ 171 PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated 172 173 Not collective 174 175 Input Parameter: 176 . q - The PetscQuadrature object 177 178 Output Parameter: 179 . Nc - The number of components 180 181 Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components. 182 183 Level: intermediate 184 185 .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData() 186 @*/ 187 PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc) 188 { 189 PetscFunctionBegin; 190 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 191 PetscValidPointer(Nc, 2); 192 *Nc = q->Nc; 193 PetscFunctionReturn(0); 194 } 195 196 /*@ 197 PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated 198 199 Not collective 200 201 Input Parameters: 202 + q - The PetscQuadrature object 203 - Nc - The number of components 204 205 Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components. 206 207 Level: intermediate 208 209 .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData() 210 @*/ 211 PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc) 212 { 213 PetscFunctionBegin; 214 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 215 q->Nc = Nc; 216 PetscFunctionReturn(0); 217 } 218 219 /*@C 220 PetscQuadratureGetData - Returns the data defining the quadrature 221 222 Not collective 223 224 Input Parameter: 225 . q - The PetscQuadrature object 226 227 Output Parameters: 228 + dim - The spatial dimension 229 . Nc - The number of components 230 . npoints - The number of quadrature points 231 . points - The coordinates of each quadrature point 232 - weights - The weight of each quadrature point 233 234 Level: intermediate 235 236 Fortran Notes: 237 From Fortran you must call PetscQuadratureRestoreData() when you are done with the data 238 239 .seealso: PetscQuadratureCreate(), PetscQuadratureSetData() 240 @*/ 241 PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[]) 242 { 243 PetscFunctionBegin; 244 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 245 if (dim) { 246 PetscValidPointer(dim, 2); 247 *dim = q->dim; 248 } 249 if (Nc) { 250 PetscValidPointer(Nc, 3); 251 *Nc = q->Nc; 252 } 253 if (npoints) { 254 PetscValidPointer(npoints, 4); 255 *npoints = q->numPoints; 256 } 257 if (points) { 258 PetscValidPointer(points, 5); 259 *points = q->points; 260 } 261 if (weights) { 262 PetscValidPointer(weights, 6); 263 *weights = q->weights; 264 } 265 PetscFunctionReturn(0); 266 } 267 268 static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[]) 269 { 270 PetscScalar *Js, *Jinvs; 271 PetscInt i, j, k; 272 PetscBLASInt bm, bn, info; 273 PetscErrorCode ierr; 274 275 PetscFunctionBegin; 276 ierr = PetscBLASIntCast(m, &bm);CHKERRQ(ierr); 277 ierr = PetscBLASIntCast(n, &bn);CHKERRQ(ierr); 278 #if defined(PETSC_USE_COMPLEX) 279 ierr = PetscMalloc2(m*n, &Js, m*n, &Jinvs);CHKERRQ(ierr); 280 for (i = 0; i < m*n; i++) Js[i] = J[i]; 281 #else 282 Js = (PetscReal *) J; 283 Jinvs = Jinv; 284 #endif 285 if (m == n) { 286 PetscBLASInt *pivots; 287 PetscScalar *W; 288 289 ierr = PetscMalloc2(m, &pivots, m, &W);CHKERRQ(ierr); 290 291 ierr = PetscArraycpy(Jinvs, Js, m * m);CHKERRQ(ierr); 292 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info)); 293 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info); 294 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &info)); 295 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info); 296 ierr = PetscFree2(pivots, W);CHKERRQ(ierr); 297 } else if (m < n) { 298 PetscScalar *JJT; 299 PetscBLASInt *pivots; 300 PetscScalar *W; 301 302 ierr = PetscMalloc1(m*m, &JJT);CHKERRQ(ierr); 303 ierr = PetscMalloc2(m, &pivots, m, &W);CHKERRQ(ierr); 304 for (i = 0; i < m; i++) { 305 for (j = 0; j < m; j++) { 306 PetscScalar val = 0.; 307 308 for (k = 0; k < n; k++) val += Js[i * n + k] * Js[j * n + k]; 309 JJT[i * m + j] = val; 310 } 311 } 312 313 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info)); 314 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info); 315 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info)); 316 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info); 317 for (i = 0; i < n; i++) { 318 for (j = 0; j < m; j++) { 319 PetscScalar val = 0.; 320 321 for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j]; 322 Jinvs[i * m + j] = val; 323 } 324 } 325 ierr = PetscFree2(pivots, W);CHKERRQ(ierr); 326 ierr = PetscFree(JJT);CHKERRQ(ierr); 327 } else { 328 PetscScalar *JTJ; 329 PetscBLASInt *pivots; 330 PetscScalar *W; 331 332 ierr = PetscMalloc1(n*n, &JTJ);CHKERRQ(ierr); 333 ierr = PetscMalloc2(n, &pivots, n, &W);CHKERRQ(ierr); 334 for (i = 0; i < n; i++) { 335 for (j = 0; j < n; j++) { 336 PetscScalar val = 0.; 337 338 for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j]; 339 JTJ[i * n + j] = val; 340 } 341 } 342 343 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bm, pivots, &info)); 344 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info); 345 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info)); 346 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info); 347 for (i = 0; i < n; i++) { 348 for (j = 0; j < m; j++) { 349 PetscScalar val = 0.; 350 351 for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k]; 352 Jinvs[i * m + j] = val; 353 } 354 } 355 ierr = PetscFree2(pivots, W);CHKERRQ(ierr); 356 ierr = PetscFree(JTJ);CHKERRQ(ierr); 357 } 358 #if defined(PETSC_USE_COMPLEX) 359 for (i = 0; i < m*n; i++) Jinv[i] = PetscRealPart(Jinvs[i]); 360 ierr = PetscFree2(Js, Jinvs);CHKERRQ(ierr); 361 #endif 362 PetscFunctionReturn(0); 363 } 364 365 /*@ 366 PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation. 367 368 Collecive on PetscQuadrature 369 370 Input Arguments: 371 + q - the quadrature functional 372 . imageDim - the dimension of the image of the transformation 373 . origin - a point in the original space 374 . originImage - the image of the origin under the transformation 375 . J - the Jacobian of the image: an [imageDim x dim] matrix in row major order 376 - 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] 377 378 Output Arguments: 379 . 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. 380 381 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. 382 383 .seealso: PetscDTAltVPullback(), PetscDTAltVPullbackMatrix() 384 @*/ 385 PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq) 386 { 387 PetscInt dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c; 388 const PetscReal *points; 389 const PetscReal *weights; 390 PetscReal *imagePoints, *imageWeights; 391 PetscReal *Jinv; 392 PetscReal *Jinvstar; 393 PetscErrorCode ierr; 394 395 PetscFunctionBegin; 396 PetscValidHeaderSpecific(q, PETSC_OBJECT_CLASSID, 1); 397 if (imageDim < PetscAbsInt(formDegree)) SETERRQ2(PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Cannot represent a %D-form in %D dimensions", PetscAbsInt(formDegree), imageDim); 398 ierr = PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights);CHKERRQ(ierr); 399 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize);CHKERRQ(ierr); 400 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); 401 Ncopies = Nc / formSize; 402 ierr = PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize);CHKERRQ(ierr); 403 imageNc = Ncopies * imageFormSize; 404 ierr = PetscMalloc1(Npoints * imageDim, &imagePoints);CHKERRQ(ierr); 405 ierr = PetscMalloc1(Npoints * imageNc, &imageWeights);CHKERRQ(ierr); 406 ierr = PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar);CHKERRQ(ierr); 407 ierr = PetscDTJacobianInverse_Internal(dim, imageDim, J, Jinv);CHKERRQ(ierr); 408 ierr = PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar);CHKERRQ(ierr); 409 for (pt = 0; pt < Npoints; pt++) { 410 const PetscReal *point = &points[pt * dim]; 411 PetscReal *imagePoint = &imagePoints[pt * imageDim]; 412 413 for (i = 0; i < imageDim; i++) { 414 PetscReal val = originImage[i]; 415 416 for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]); 417 imagePoint[i] = val; 418 } 419 for (c = 0; c < Ncopies; c++) { 420 const PetscReal *form = &weights[pt * Nc + c * formSize]; 421 PetscReal *imageForm = &imageWeights[pt * imageNc + c * imageFormSize]; 422 423 for (i = 0; i < imageFormSize; i++) { 424 PetscReal val = 0.; 425 426 for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j]; 427 imageForm[i] = val; 428 } 429 } 430 } 431 ierr = PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq);CHKERRQ(ierr); 432 ierr = PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights);CHKERRQ(ierr); 433 ierr = PetscFree2(Jinv, Jinvstar);CHKERRQ(ierr); 434 PetscFunctionReturn(0); 435 } 436 437 /*@C 438 PetscQuadratureSetData - Sets the data defining the quadrature 439 440 Not collective 441 442 Input Parameters: 443 + q - The PetscQuadrature object 444 . dim - The spatial dimension 445 . Nc - The number of components 446 . npoints - The number of quadrature points 447 . points - The coordinates of each quadrature point 448 - weights - The weight of each quadrature point 449 450 Note: This routine owns the references to points and weights, so they must be allocated using PetscMalloc() and the user should not free them. 451 452 Level: intermediate 453 454 .seealso: PetscQuadratureCreate(), PetscQuadratureGetData() 455 @*/ 456 PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[]) 457 { 458 PetscFunctionBegin; 459 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 460 if (dim >= 0) q->dim = dim; 461 if (Nc >= 0) q->Nc = Nc; 462 if (npoints >= 0) q->numPoints = npoints; 463 if (points) { 464 PetscValidPointer(points, 4); 465 q->points = points; 466 } 467 if (weights) { 468 PetscValidPointer(weights, 5); 469 q->weights = weights; 470 } 471 PetscFunctionReturn(0); 472 } 473 474 static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v) 475 { 476 PetscInt q, d, c; 477 PetscViewerFormat format; 478 PetscErrorCode ierr; 479 480 PetscFunctionBegin; 481 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);} 482 else {ierr = PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D)\n", quad->order, quad->numPoints, quad->dim);CHKERRQ(ierr);} 483 ierr = PetscViewerGetFormat(v, &format);CHKERRQ(ierr); 484 if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 485 for (q = 0; q < quad->numPoints; ++q) { 486 ierr = PetscViewerASCIIPrintf(v, "p%D (", q);CHKERRQ(ierr); 487 ierr = PetscViewerASCIIUseTabs(v, PETSC_FALSE);CHKERRQ(ierr); 488 for (d = 0; d < quad->dim; ++d) { 489 if (d) {ierr = PetscViewerASCIIPrintf(v, ", ");CHKERRQ(ierr);} 490 ierr = PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q*quad->dim+d]);CHKERRQ(ierr); 491 } 492 ierr = PetscViewerASCIIPrintf(v, ") ");CHKERRQ(ierr); 493 if (quad->Nc > 1) {ierr = PetscViewerASCIIPrintf(v, "w%D (", q);CHKERRQ(ierr);} 494 for (c = 0; c < quad->Nc; ++c) { 495 if (c) {ierr = PetscViewerASCIIPrintf(v, ", ");CHKERRQ(ierr);} 496 ierr = PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q*quad->Nc+c]);CHKERRQ(ierr); 497 } 498 if (quad->Nc > 1) {ierr = PetscViewerASCIIPrintf(v, ")");CHKERRQ(ierr);} 499 ierr = PetscViewerASCIIPrintf(v, "\n");CHKERRQ(ierr); 500 ierr = PetscViewerASCIIUseTabs(v, PETSC_TRUE);CHKERRQ(ierr); 501 } 502 PetscFunctionReturn(0); 503 } 504 505 /*@C 506 PetscQuadratureView - Views a PetscQuadrature object 507 508 Collective on quad 509 510 Input Parameters: 511 + quad - The PetscQuadrature object 512 - viewer - The PetscViewer object 513 514 Level: beginner 515 516 .seealso: PetscQuadratureCreate(), PetscQuadratureGetData() 517 @*/ 518 PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer) 519 { 520 PetscBool iascii; 521 PetscErrorCode ierr; 522 523 PetscFunctionBegin; 524 PetscValidHeader(quad, 1); 525 if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 526 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) quad), &viewer);CHKERRQ(ierr);} 527 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 528 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 529 if (iascii) {ierr = PetscQuadratureView_Ascii(quad, viewer);CHKERRQ(ierr);} 530 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 531 PetscFunctionReturn(0); 532 } 533 534 /*@C 535 PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement 536 537 Not collective 538 539 Input Parameter: 540 + q - The original PetscQuadrature 541 . numSubelements - The number of subelements the original element is divided into 542 . v0 - An array of the initial points for each subelement 543 - jac - An array of the Jacobian mappings from the reference to each subelement 544 545 Output Parameters: 546 . dim - The dimension 547 548 Note: Together v0 and jac define an affine mapping from the original reference element to each subelement 549 550 Not available from Fortran 551 552 Level: intermediate 553 554 .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension() 555 @*/ 556 PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref) 557 { 558 const PetscReal *points, *weights; 559 PetscReal *pointsRef, *weightsRef; 560 PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e; 561 PetscErrorCode ierr; 562 563 PetscFunctionBegin; 564 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1); 565 PetscValidPointer(v0, 3); 566 PetscValidPointer(jac, 4); 567 PetscValidPointer(qref, 5); 568 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, qref);CHKERRQ(ierr); 569 ierr = PetscQuadratureGetOrder(q, &order);CHKERRQ(ierr); 570 ierr = PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);CHKERRQ(ierr); 571 npointsRef = npoints*numSubelements; 572 ierr = PetscMalloc1(npointsRef*dim,&pointsRef);CHKERRQ(ierr); 573 ierr = PetscMalloc1(npointsRef*Nc, &weightsRef);CHKERRQ(ierr); 574 for (c = 0; c < numSubelements; ++c) { 575 for (p = 0; p < npoints; ++p) { 576 for (d = 0; d < dim; ++d) { 577 pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d]; 578 for (e = 0; e < dim; ++e) { 579 pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0); 580 } 581 } 582 /* Could also use detJ here */ 583 for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements; 584 } 585 } 586 ierr = PetscQuadratureSetOrder(*qref, order);CHKERRQ(ierr); 587 ierr = PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);CHKERRQ(ierr); 588 PetscFunctionReturn(0); 589 } 590 591 /*@ 592 PetscDTLegendreEval - evaluate Legendre polynomial at points 593 594 Not Collective 595 596 Input Arguments: 597 + npoints - number of spatial points to evaluate at 598 . points - array of locations to evaluate at 599 . ndegree - number of basis degrees to evaluate 600 - degrees - sorted array of degrees to evaluate 601 602 Output Arguments: 603 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL) 604 . D - row-oriented derivative evaluation matrix (or NULL) 605 - D2 - row-oriented second derivative evaluation matrix (or NULL) 606 607 Level: intermediate 608 609 .seealso: PetscDTGaussQuadrature() 610 @*/ 611 PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2) 612 { 613 PetscInt i,maxdegree; 614 615 PetscFunctionBegin; 616 if (!npoints || !ndegree) PetscFunctionReturn(0); 617 maxdegree = degrees[ndegree-1]; 618 for (i=0; i<npoints; i++) { 619 PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x; 620 PetscInt j,k; 621 x = points[i]; 622 pm2 = 0; 623 pm1 = 1; 624 pd2 = 0; 625 pd1 = 0; 626 pdd2 = 0; 627 pdd1 = 0; 628 k = 0; 629 if (degrees[k] == 0) { 630 if (B) B[i*ndegree+k] = pm1; 631 if (D) D[i*ndegree+k] = pd1; 632 if (D2) D2[i*ndegree+k] = pdd1; 633 k++; 634 } 635 for (j=1; j<=maxdegree; j++,k++) { 636 PetscReal p,d,dd; 637 p = ((2*j-1)*x*pm1 - (j-1)*pm2)/j; 638 d = pd2 + (2*j-1)*pm1; 639 dd = pdd2 + (2*j-1)*pd1; 640 pm2 = pm1; 641 pm1 = p; 642 pd2 = pd1; 643 pd1 = d; 644 pdd2 = pdd1; 645 pdd1 = dd; 646 if (degrees[k] == j) { 647 if (B) B[i*ndegree+k] = p; 648 if (D) D[i*ndegree+k] = d; 649 if (D2) D2[i*ndegree+k] = dd; 650 } 651 } 652 } 653 PetscFunctionReturn(0); 654 } 655 656 /*@ 657 PetscDTGaussQuadrature - create Gauss quadrature 658 659 Not Collective 660 661 Input Arguments: 662 + npoints - number of points 663 . a - left end of interval (often-1) 664 - b - right end of interval (often +1) 665 666 Output Arguments: 667 + x - quadrature points 668 - w - quadrature weights 669 670 Level: intermediate 671 672 References: 673 . 1. - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969. 674 675 .seealso: PetscDTLegendreEval() 676 @*/ 677 PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w) 678 { 679 PetscErrorCode ierr; 680 PetscInt i; 681 PetscReal *work; 682 PetscScalar *Z; 683 PetscBLASInt N,LDZ,info; 684 685 PetscFunctionBegin; 686 ierr = PetscCitationsRegister(GaussCitation, &GaussCite);CHKERRQ(ierr); 687 /* Set up the Golub-Welsch system */ 688 for (i=0; i<npoints; i++) { 689 x[i] = 0; /* diagonal is 0 */ 690 if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i)); 691 } 692 ierr = PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);CHKERRQ(ierr); 693 ierr = PetscBLASIntCast(npoints,&N);CHKERRQ(ierr); 694 LDZ = N; 695 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 696 #if defined(PETSC_MISSING_LAPACK_STEQR) 697 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"STEQR - Lapack routine is unavailable."); 698 #else 699 PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info)); 700 #endif 701 ierr = PetscFPTrapPop();CHKERRQ(ierr); 702 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error"); 703 704 for (i=0; i<(npoints+1)/2; i++) { 705 PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */ 706 x[i] = (a+b)/2 - y*(b-a)/2; 707 if (x[i] == -0.0) x[i] = 0.0; 708 x[npoints-i-1] = (a+b)/2 + y*(b-a)/2; 709 710 w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints]))); 711 } 712 ierr = PetscFree2(Z,work);CHKERRQ(ierr); 713 PetscFunctionReturn(0); 714 } 715 716 static void qAndLEvaluation(PetscInt n, PetscReal x, PetscReal *q, PetscReal *qp, PetscReal *Ln) 717 /* 718 Compute the polynomial q(x) = L_{N+1}(x) - L_{n-1}(x) and its derivative in 719 addition to L_N(x) as these are needed for computing the GLL points via Newton's method. 720 Reference: "Implementing Spectral Methods for Partial Differential Equations: Algorithms 721 for Scientists and Engineers" by David A. Kopriva. 722 */ 723 { 724 PetscInt k; 725 726 PetscReal Lnp; 727 PetscReal Lnp1, Lnp1p; 728 PetscReal Lnm1, Lnm1p; 729 PetscReal Lnm2, Lnm2p; 730 731 Lnm1 = 1.0; 732 *Ln = x; 733 Lnm1p = 0.0; 734 Lnp = 1.0; 735 736 for (k=2; k<=n; ++k) { 737 Lnm2 = Lnm1; 738 Lnm1 = *Ln; 739 Lnm2p = Lnm1p; 740 Lnm1p = Lnp; 741 *Ln = (2.*((PetscReal)k)-1.)/(1.0*((PetscReal)k))*x*Lnm1 - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm2; 742 Lnp = Lnm2p + (2.0*((PetscReal)k)-1.)*Lnm1; 743 } 744 k = n+1; 745 Lnp1 = (2.*((PetscReal)k)-1.)/(((PetscReal)k))*x*(*Ln) - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm1; 746 Lnp1p = Lnm1p + (2.0*((PetscReal)k)-1.)*(*Ln); 747 *q = Lnp1 - Lnm1; 748 *qp = Lnp1p - Lnm1p; 749 } 750 751 /*@C 752 PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre 753 nodes of a given size on the domain [-1,1] 754 755 Not Collective 756 757 Input Parameter: 758 + n - number of grid nodes 759 - type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON 760 761 Output Arguments: 762 + x - quadrature points 763 - w - quadrature weights 764 765 Notes: 766 For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not 767 close enough to the desired solution 768 769 These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes 770 771 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 772 773 Level: intermediate 774 775 .seealso: PetscDTGaussQuadrature() 776 777 @*/ 778 PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal *x,PetscReal *w) 779 { 780 PetscErrorCode ierr; 781 782 PetscFunctionBegin; 783 if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide at least 2 grid points per element"); 784 785 if (type == PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA) { 786 PetscReal *M,si; 787 PetscBLASInt bn,lierr; 788 PetscReal x0,z0,z1,z2; 789 PetscInt i,p = npoints - 1,nn; 790 791 x[0] =-1.0; 792 x[npoints-1] = 1.0; 793 if (npoints-2 > 0){ 794 ierr = PetscMalloc1(npoints-1,&M);CHKERRQ(ierr); 795 for (i=0; i<npoints-2; i++) { 796 si = ((PetscReal)i)+1.0; 797 M[i]=0.5*PetscSqrtReal(si*(si+2.0)/((si+0.5)*(si+1.5))); 798 } 799 ierr = PetscBLASIntCast(npoints-2,&bn);CHKERRQ(ierr); 800 ierr = PetscArrayzero(&x[1],bn);CHKERRQ(ierr); 801 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 802 x0=0; 803 #if defined(PETSC_MISSING_LAPACK_STEQR) 804 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"STEQR - Lapack routine is unavailable."); 805 #else 806 PetscStackCallBLAS("LAPACKsteqr",LAPACKREALsteqr_("N",&bn,&x[1],M,&x0,&bn,M,&lierr)); 807 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in STERF Lapack routine %d",(int)lierr); 808 #endif 809 ierr = PetscFPTrapPop();CHKERRQ(ierr); 810 ierr = PetscFree(M);CHKERRQ(ierr); 811 } 812 if ((npoints-1)%2==0) { 813 x[(npoints-1)/2] = 0.0; /* hard wire to exactly 0.0 since linear algebra produces nonzero */ 814 } 815 816 w[0] = w[p] = 2.0/(((PetscReal)(p))*(((PetscReal)p)+1.0)); 817 z2 = -1.; /* Dummy value to avoid -Wmaybe-initialized */ 818 for (i=1; i<p; i++) { 819 x0 = x[i]; 820 z0 = 1.0; 821 z1 = x0; 822 for (nn=1; nn<p; nn++) { 823 z2 = x0*z1*(2.0*((PetscReal)nn)+1.0)/(((PetscReal)nn)+1.0)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.0)); 824 z0 = z1; 825 z1 = z2; 826 } 827 w[i]=2.0/(((PetscReal)p)*(((PetscReal)p)+1.0)*z2*z2); 828 } 829 } else { 830 PetscInt j,m; 831 PetscReal z1,z,q,qp,Ln; 832 PetscReal *pt; 833 ierr = PetscMalloc1(npoints,&pt);CHKERRQ(ierr); 834 835 if (npoints > 30) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON produces incorrect answers for n > 30"); 836 x[0] = -1.0; 837 x[npoints-1] = 1.0; 838 w[0] = w[npoints-1] = 2./(((PetscReal)npoints)*(((PetscReal)npoints)-1.0)); 839 m = (npoints-1)/2; /* The roots are symmetric, so we only find half of them. */ 840 for (j=1; j<=m; j++) { /* Loop over the desired roots. */ 841 z = -1.0*PetscCosReal((PETSC_PI*((PetscReal)j)+0.25)/(((PetscReal)npoints)-1.0))-(3.0/(8.0*(((PetscReal)npoints)-1.0)*PETSC_PI))*(1.0/(((PetscReal)j)+0.25)); 842 /* Starting with the above approximation to the ith root, we enter */ 843 /* the main loop of refinement by Newton's method. */ 844 do { 845 qAndLEvaluation(npoints-1,z,&q,&qp,&Ln); 846 z1 = z; 847 z = z1-q/qp; /* Newton's method. */ 848 } while (PetscAbs(z-z1) > 10.*PETSC_MACHINE_EPSILON); 849 qAndLEvaluation(npoints-1,z,&q,&qp,&Ln); 850 851 x[j] = z; 852 x[npoints-1-j] = -z; /* and put in its symmetric counterpart. */ 853 w[j] = 2.0/(((PetscReal)npoints)*(((PetscReal)npoints)-1.)*Ln*Ln); /* Compute the weight */ 854 w[npoints-1-j] = w[j]; /* and its symmetric counterpart. */ 855 pt[j]=qp; 856 } 857 858 if ((npoints-1)%2==0) { 859 qAndLEvaluation(npoints-1,0.0,&q,&qp,&Ln); 860 x[(npoints-1)/2] = 0.0; 861 w[(npoints-1)/2] = 2.0/(((PetscReal)npoints)*(((PetscReal)npoints)-1.)*Ln*Ln); 862 } 863 ierr = PetscFree(pt);CHKERRQ(ierr); 864 } 865 PetscFunctionReturn(0); 866 } 867 868 /*@ 869 PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature 870 871 Not Collective 872 873 Input Arguments: 874 + dim - The spatial dimension 875 . Nc - The number of components 876 . npoints - number of points in one dimension 877 . a - left end of interval (often-1) 878 - b - right end of interval (often +1) 879 880 Output Argument: 881 . q - A PetscQuadrature object 882 883 Level: intermediate 884 885 .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval() 886 @*/ 887 PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 888 { 889 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c; 890 PetscReal *x, *w, *xw, *ww; 891 PetscErrorCode ierr; 892 893 PetscFunctionBegin; 894 ierr = PetscMalloc1(totpoints*dim,&x);CHKERRQ(ierr); 895 ierr = PetscMalloc1(totpoints*Nc,&w);CHKERRQ(ierr); 896 /* Set up the Golub-Welsch system */ 897 switch (dim) { 898 case 0: 899 ierr = PetscFree(x);CHKERRQ(ierr); 900 ierr = PetscFree(w);CHKERRQ(ierr); 901 ierr = PetscMalloc1(1, &x);CHKERRQ(ierr); 902 ierr = PetscMalloc1(Nc, &w);CHKERRQ(ierr); 903 x[0] = 0.0; 904 for (c = 0; c < Nc; ++c) w[c] = 1.0; 905 break; 906 case 1: 907 ierr = PetscMalloc1(npoints,&ww);CHKERRQ(ierr); 908 ierr = PetscDTGaussQuadrature(npoints, a, b, x, ww);CHKERRQ(ierr); 909 for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i]; 910 ierr = PetscFree(ww);CHKERRQ(ierr); 911 break; 912 case 2: 913 ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr); 914 ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr); 915 for (i = 0; i < npoints; ++i) { 916 for (j = 0; j < npoints; ++j) { 917 x[(i*npoints+j)*dim+0] = xw[i]; 918 x[(i*npoints+j)*dim+1] = xw[j]; 919 for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j]; 920 } 921 } 922 ierr = PetscFree2(xw,ww);CHKERRQ(ierr); 923 break; 924 case 3: 925 ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr); 926 ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr); 927 for (i = 0; i < npoints; ++i) { 928 for (j = 0; j < npoints; ++j) { 929 for (k = 0; k < npoints; ++k) { 930 x[((i*npoints+j)*npoints+k)*dim+0] = xw[i]; 931 x[((i*npoints+j)*npoints+k)*dim+1] = xw[j]; 932 x[((i*npoints+j)*npoints+k)*dim+2] = xw[k]; 933 for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k]; 934 } 935 } 936 } 937 ierr = PetscFree2(xw,ww);CHKERRQ(ierr); 938 break; 939 default: 940 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim); 941 } 942 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 943 ierr = PetscQuadratureSetOrder(*q, 2*npoints-1);CHKERRQ(ierr); 944 ierr = PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);CHKERRQ(ierr); 945 ierr = PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");CHKERRQ(ierr); 946 PetscFunctionReturn(0); 947 } 948 949 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x. 950 Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */ 951 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 952 { 953 PetscReal apb, pn1, pn2; 954 PetscInt k; 955 956 PetscFunctionBegin; 957 if (!n) {*P = 1.0; PetscFunctionReturn(0);} 958 if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); PetscFunctionReturn(0);} 959 apb = a + b; 960 pn2 = 1.0; 961 pn1 = 0.5 * (a - b + (apb + 2.0) * x); 962 *P = 0.0; 963 for (k = 2; k < n+1; ++k) { 964 PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0); 965 PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b); 966 PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb); 967 PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb); 968 969 a2 = a2 / a1; 970 a3 = a3 / a1; 971 a4 = a4 / a1; 972 *P = (a2 + a3 * x) * pn1 - a4 * pn2; 973 pn2 = pn1; 974 pn1 = *P; 975 } 976 PetscFunctionReturn(0); 977 } 978 979 /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */ 980 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 981 { 982 PetscReal nP; 983 PetscErrorCode ierr; 984 985 PetscFunctionBegin; 986 if (!n) {*P = 0.0; PetscFunctionReturn(0);} 987 ierr = PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);CHKERRQ(ierr); 988 *P = 0.5 * (a + b + n + 1) * nP; 989 PetscFunctionReturn(0); 990 } 991 992 /* Maps from [-1,1]^2 to the (-1,1) reference triangle */ 993 PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta) 994 { 995 PetscFunctionBegin; 996 *xi = 0.5 * (1.0 + x) * (1.0 - y) - 1.0; 997 *eta = y; 998 PetscFunctionReturn(0); 999 } 1000 1001 /* Maps from [-1,1]^2 to the (-1,1) reference triangle */ 1002 PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta) 1003 { 1004 PetscFunctionBegin; 1005 *xi = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0; 1006 *eta = 0.5 * (1.0 + y) * (1.0 - z) - 1.0; 1007 *zeta = z; 1008 PetscFunctionReturn(0); 1009 } 1010 1011 static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w) 1012 { 1013 PetscInt maxIter = 100; 1014 PetscReal eps = 1.0e-8; 1015 PetscReal a1, a2, a3, a4, a5, a6; 1016 PetscInt k; 1017 PetscErrorCode ierr; 1018 1019 PetscFunctionBegin; 1020 1021 a1 = PetscPowReal(2.0, a+b+1); 1022 #if defined(PETSC_HAVE_TGAMMA) 1023 a2 = PetscTGamma(a + npoints + 1); 1024 a3 = PetscTGamma(b + npoints + 1); 1025 a4 = PetscTGamma(a + b + npoints + 1); 1026 #else 1027 { 1028 PetscInt ia, ib; 1029 1030 ia = (PetscInt) a; 1031 ib = (PetscInt) b; 1032 if (ia == a && ib == b && ia + npoints + 1 > 0 && ib + npoints + 1 > 0 && ia + ib + npoints + 1 > 0) { /* All gamma(x) terms are (x-1)! terms */ 1033 ierr = PetscDTFactorial(ia + npoints, &a2);CHKERRQ(ierr); 1034 ierr = PetscDTFactorial(ib + npoints, &a3);CHKERRQ(ierr); 1035 ierr = PetscDTFactorial(ia + ib + npoints, &a4);CHKERRQ(ierr); 1036 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable."); 1037 } 1038 #endif 1039 1040 ierr = PetscDTFactorial(npoints, &a5);CHKERRQ(ierr); 1041 a6 = a1 * a2 * a3 / a4 / a5; 1042 /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses. 1043 Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */ 1044 for (k = 0; k < npoints; ++k) { 1045 PetscReal r = -PetscCosReal((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP; 1046 PetscInt j; 1047 1048 if (k > 0) r = 0.5 * (r + x[k-1]); 1049 for (j = 0; j < maxIter; ++j) { 1050 PetscReal s = 0.0, delta, f, fp; 1051 PetscInt i; 1052 1053 for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]); 1054 ierr = PetscDTComputeJacobi(a, b, npoints, r, &f);CHKERRQ(ierr); 1055 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);CHKERRQ(ierr); 1056 delta = f / (fp - f * s); 1057 r = r - delta; 1058 if (PetscAbsReal(delta) < eps) break; 1059 } 1060 x[k] = r; 1061 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);CHKERRQ(ierr); 1062 w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP); 1063 } 1064 PetscFunctionReturn(0); 1065 } 1066 1067 /*@ 1068 PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex 1069 1070 Not Collective 1071 1072 Input Arguments: 1073 + dim - The simplex dimension 1074 . Nc - The number of components 1075 . npoints - The number of points in one dimension 1076 . a - left end of interval (often-1) 1077 - b - right end of interval (often +1) 1078 1079 Output Argument: 1080 . q - A PetscQuadrature object 1081 1082 Level: intermediate 1083 1084 References: 1085 . 1. - Karniadakis and Sherwin. FIAT 1086 1087 .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature() 1088 @*/ 1089 PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 1090 { 1091 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints; 1092 PetscReal *px, *wx, *py, *wy, *pz, *wz, *x, *w; 1093 PetscInt i, j, k, c; 1094 PetscErrorCode ierr; 1095 1096 PetscFunctionBegin; 1097 if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now"); 1098 ierr = PetscMalloc1(totpoints*dim, &x);CHKERRQ(ierr); 1099 ierr = PetscMalloc1(totpoints*Nc, &w);CHKERRQ(ierr); 1100 switch (dim) { 1101 case 0: 1102 ierr = PetscFree(x);CHKERRQ(ierr); 1103 ierr = PetscFree(w);CHKERRQ(ierr); 1104 ierr = PetscMalloc1(1, &x);CHKERRQ(ierr); 1105 ierr = PetscMalloc1(Nc, &w);CHKERRQ(ierr); 1106 x[0] = 0.0; 1107 for (c = 0; c < Nc; ++c) w[c] = 1.0; 1108 break; 1109 case 1: 1110 ierr = PetscMalloc1(npoints,&wx);CHKERRQ(ierr); 1111 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, wx);CHKERRQ(ierr); 1112 for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i]; 1113 ierr = PetscFree(wx);CHKERRQ(ierr); 1114 break; 1115 case 2: 1116 ierr = PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);CHKERRQ(ierr); 1117 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);CHKERRQ(ierr); 1118 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);CHKERRQ(ierr); 1119 for (i = 0; i < npoints; ++i) { 1120 for (j = 0; j < npoints; ++j) { 1121 ierr = PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);CHKERRQ(ierr); 1122 for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j]; 1123 } 1124 } 1125 ierr = PetscFree4(px,wx,py,wy);CHKERRQ(ierr); 1126 break; 1127 case 3: 1128 ierr = PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);CHKERRQ(ierr); 1129 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);CHKERRQ(ierr); 1130 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);CHKERRQ(ierr); 1131 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);CHKERRQ(ierr); 1132 for (i = 0; i < npoints; ++i) { 1133 for (j = 0; j < npoints; ++j) { 1134 for (k = 0; k < npoints; ++k) { 1135 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); 1136 for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k]; 1137 } 1138 } 1139 } 1140 ierr = PetscFree6(px,wx,py,wy,pz,wz);CHKERRQ(ierr); 1141 break; 1142 default: 1143 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim); 1144 } 1145 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 1146 ierr = PetscQuadratureSetOrder(*q, 2*npoints-1);CHKERRQ(ierr); 1147 ierr = PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);CHKERRQ(ierr); 1148 ierr = PetscObjectChangeTypeName((PetscObject)*q,"GaussJacobi");CHKERRQ(ierr); 1149 PetscFunctionReturn(0); 1150 } 1151 1152 /*@ 1153 PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell 1154 1155 Not Collective 1156 1157 Input Arguments: 1158 + dim - The cell dimension 1159 . level - The number of points in one dimension, 2^l 1160 . a - left end of interval (often-1) 1161 - b - right end of interval (often +1) 1162 1163 Output Argument: 1164 . q - A PetscQuadrature object 1165 1166 Level: intermediate 1167 1168 .seealso: PetscDTGaussTensorQuadrature() 1169 @*/ 1170 PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q) 1171 { 1172 const PetscInt p = 16; /* Digits of precision in the evaluation */ 1173 const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */ 1174 const PetscReal beta = (b+a)/2.; /* Center of the integration interval */ 1175 const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */ 1176 PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */ 1177 PetscReal wk = 0.5*PETSC_PI; /* Quadrature weight at x_k */ 1178 PetscReal *x, *w; 1179 PetscInt K, k, npoints; 1180 PetscErrorCode ierr; 1181 1182 PetscFunctionBegin; 1183 if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim); 1184 if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits"); 1185 /* Find K such that the weights are < 32 digits of precision */ 1186 for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) { 1187 wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h))); 1188 } 1189 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 1190 ierr = PetscQuadratureSetOrder(*q, 2*K+1);CHKERRQ(ierr); 1191 npoints = 2*K-1; 1192 ierr = PetscMalloc1(npoints*dim, &x);CHKERRQ(ierr); 1193 ierr = PetscMalloc1(npoints, &w);CHKERRQ(ierr); 1194 /* Center term */ 1195 x[0] = beta; 1196 w[0] = 0.5*alpha*PETSC_PI; 1197 for (k = 1; k < K; ++k) { 1198 wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1199 xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h)); 1200 x[2*k-1] = -alpha*xk+beta; 1201 w[2*k-1] = wk; 1202 x[2*k+0] = alpha*xk+beta; 1203 w[2*k+0] = wk; 1204 } 1205 ierr = PetscQuadratureSetData(*q, dim, 1, npoints, x, w);CHKERRQ(ierr); 1206 PetscFunctionReturn(0); 1207 } 1208 1209 PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 1210 { 1211 const PetscInt p = 16; /* Digits of precision in the evaluation */ 1212 const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */ 1213 const PetscReal beta = (b+a)/2.; /* Center of the integration interval */ 1214 PetscReal h = 1.0; /* Step size, length between x_k */ 1215 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 1216 PetscReal osum = 0.0; /* Integral on last level */ 1217 PetscReal psum = 0.0; /* Integral on the level before the last level */ 1218 PetscReal sum; /* Integral on current level */ 1219 PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 1220 PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 1221 PetscReal wk; /* Quadrature weight at x_k */ 1222 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 1223 PetscInt d; /* Digits of precision in the integral */ 1224 1225 PetscFunctionBegin; 1226 if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 1227 /* Center term */ 1228 func(beta, &lval); 1229 sum = 0.5*alpha*PETSC_PI*lval; 1230 /* */ 1231 do { 1232 PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4; 1233 PetscInt k = 1; 1234 1235 ++l; 1236 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */ 1237 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 1238 psum = osum; 1239 osum = sum; 1240 h *= 0.5; 1241 sum *= 0.5; 1242 do { 1243 wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1244 yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1245 lx = -alpha*(1.0 - yk)+beta; 1246 rx = alpha*(1.0 - yk)+beta; 1247 func(lx, &lval); 1248 func(rx, &rval); 1249 lterm = alpha*wk*lval; 1250 maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm); 1251 sum += lterm; 1252 rterm = alpha*wk*rval; 1253 maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm); 1254 sum += rterm; 1255 ++k; 1256 /* Only need to evaluate every other point on refined levels */ 1257 if (l != 1) ++k; 1258 } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */ 1259 1260 d1 = PetscLog10Real(PetscAbsReal(sum - osum)); 1261 d2 = PetscLog10Real(PetscAbsReal(sum - psum)); 1262 d3 = PetscLog10Real(maxTerm) - p; 1263 if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0; 1264 else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm))); 1265 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4))); 1266 } while (d < digits && l < 12); 1267 *sol = sum; 1268 1269 PetscFunctionReturn(0); 1270 } 1271 1272 #if defined(PETSC_HAVE_MPFR) 1273 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 1274 { 1275 const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */ 1276 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 1277 mpfr_t alpha; /* Half-width of the integration interval */ 1278 mpfr_t beta; /* Center of the integration interval */ 1279 mpfr_t h; /* Step size, length between x_k */ 1280 mpfr_t osum; /* Integral on last level */ 1281 mpfr_t psum; /* Integral on the level before the last level */ 1282 mpfr_t sum; /* Integral on current level */ 1283 mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 1284 mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 1285 mpfr_t wk; /* Quadrature weight at x_k */ 1286 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 1287 PetscInt d; /* Digits of precision in the integral */ 1288 mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp; 1289 1290 PetscFunctionBegin; 1291 if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 1292 /* Create high precision storage */ 1293 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); 1294 /* Initialization */ 1295 mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN); 1296 mpfr_set_d(beta, 0.5*(b+a), MPFR_RNDN); 1297 mpfr_set_d(osum, 0.0, MPFR_RNDN); 1298 mpfr_set_d(psum, 0.0, MPFR_RNDN); 1299 mpfr_set_d(h, 1.0, MPFR_RNDN); 1300 mpfr_const_pi(pi2, MPFR_RNDN); 1301 mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN); 1302 /* Center term */ 1303 func(0.5*(b+a), &lval); 1304 mpfr_set(sum, pi2, MPFR_RNDN); 1305 mpfr_mul(sum, sum, alpha, MPFR_RNDN); 1306 mpfr_mul_d(sum, sum, lval, MPFR_RNDN); 1307 /* */ 1308 do { 1309 PetscReal d1, d2, d3, d4; 1310 PetscInt k = 1; 1311 1312 ++l; 1313 mpfr_set_d(maxTerm, 0.0, MPFR_RNDN); 1314 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */ 1315 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 1316 mpfr_set(psum, osum, MPFR_RNDN); 1317 mpfr_set(osum, sum, MPFR_RNDN); 1318 mpfr_mul_d(h, h, 0.5, MPFR_RNDN); 1319 mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN); 1320 do { 1321 mpfr_set_si(kh, k, MPFR_RNDN); 1322 mpfr_mul(kh, kh, h, MPFR_RNDN); 1323 /* Weight */ 1324 mpfr_set(wk, h, MPFR_RNDN); 1325 mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN); 1326 mpfr_mul(msinh, msinh, pi2, MPFR_RNDN); 1327 mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN); 1328 mpfr_cosh(tmp, msinh, MPFR_RNDN); 1329 mpfr_sqr(tmp, tmp, MPFR_RNDN); 1330 mpfr_mul(wk, wk, mcosh, MPFR_RNDN); 1331 mpfr_div(wk, wk, tmp, MPFR_RNDN); 1332 /* Abscissa */ 1333 mpfr_set_d(yk, 1.0, MPFR_RNDZ); 1334 mpfr_cosh(tmp, msinh, MPFR_RNDN); 1335 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 1336 mpfr_exp(tmp, msinh, MPFR_RNDN); 1337 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 1338 /* Quadrature points */ 1339 mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ); 1340 mpfr_mul(lx, lx, alpha, MPFR_RNDU); 1341 mpfr_add(lx, lx, beta, MPFR_RNDU); 1342 mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ); 1343 mpfr_mul(rx, rx, alpha, MPFR_RNDD); 1344 mpfr_add(rx, rx, beta, MPFR_RNDD); 1345 /* Evaluation */ 1346 func(mpfr_get_d(lx, MPFR_RNDU), &lval); 1347 func(mpfr_get_d(rx, MPFR_RNDD), &rval); 1348 /* Update */ 1349 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 1350 mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN); 1351 mpfr_add(sum, sum, tmp, MPFR_RNDN); 1352 mpfr_abs(tmp, tmp, MPFR_RNDN); 1353 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 1354 mpfr_set(curTerm, tmp, MPFR_RNDN); 1355 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 1356 mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN); 1357 mpfr_add(sum, sum, tmp, MPFR_RNDN); 1358 mpfr_abs(tmp, tmp, MPFR_RNDN); 1359 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 1360 mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN); 1361 ++k; 1362 /* Only need to evaluate every other point on refined levels */ 1363 if (l != 1) ++k; 1364 mpfr_log10(tmp, wk, MPFR_RNDN); 1365 mpfr_abs(tmp, tmp, MPFR_RNDN); 1366 } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */ 1367 mpfr_sub(tmp, sum, osum, MPFR_RNDN); 1368 mpfr_abs(tmp, tmp, MPFR_RNDN); 1369 mpfr_log10(tmp, tmp, MPFR_RNDN); 1370 d1 = mpfr_get_d(tmp, MPFR_RNDN); 1371 mpfr_sub(tmp, sum, psum, MPFR_RNDN); 1372 mpfr_abs(tmp, tmp, MPFR_RNDN); 1373 mpfr_log10(tmp, tmp, MPFR_RNDN); 1374 d2 = mpfr_get_d(tmp, MPFR_RNDN); 1375 mpfr_log10(tmp, maxTerm, MPFR_RNDN); 1376 d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits; 1377 mpfr_log10(tmp, curTerm, MPFR_RNDN); 1378 d4 = mpfr_get_d(tmp, MPFR_RNDN); 1379 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4))); 1380 } while (d < digits && l < 8); 1381 *sol = mpfr_get_d(sum, MPFR_RNDN); 1382 /* Cleanup */ 1383 mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL); 1384 PetscFunctionReturn(0); 1385 } 1386 #else 1387 1388 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 1389 { 1390 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp"); 1391 } 1392 #endif 1393 1394 /* Overwrites A. Can only handle full-rank problems with m>=n 1395 * A in column-major format 1396 * Ainv in row-major format 1397 * tau has length m 1398 * worksize must be >= max(1,n) 1399 */ 1400 static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work) 1401 { 1402 PetscErrorCode ierr; 1403 PetscBLASInt M,N,K,lda,ldb,ldwork,info; 1404 PetscScalar *A,*Ainv,*R,*Q,Alpha; 1405 1406 PetscFunctionBegin; 1407 #if defined(PETSC_USE_COMPLEX) 1408 { 1409 PetscInt i,j; 1410 ierr = PetscMalloc2(m*n,&A,m*n,&Ainv);CHKERRQ(ierr); 1411 for (j=0; j<n; j++) { 1412 for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j]; 1413 } 1414 mstride = m; 1415 } 1416 #else 1417 A = A_in; 1418 Ainv = Ainv_out; 1419 #endif 1420 1421 ierr = PetscBLASIntCast(m,&M);CHKERRQ(ierr); 1422 ierr = PetscBLASIntCast(n,&N);CHKERRQ(ierr); 1423 ierr = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr); 1424 ierr = PetscBLASIntCast(worksize,&ldwork);CHKERRQ(ierr); 1425 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1426 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info)); 1427 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1428 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error"); 1429 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */ 1430 1431 /* Extract an explicit representation of Q */ 1432 Q = Ainv; 1433 ierr = PetscArraycpy(Q,A,mstride*n);CHKERRQ(ierr); 1434 K = N; /* full rank */ 1435 #if defined(PETSC_MISSING_LAPACK_ORGQR) 1436 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"ORGQR - Lapack routine is unavailable."); 1437 #else 1438 PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info)); 1439 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error"); 1440 #endif 1441 1442 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */ 1443 Alpha = 1.0; 1444 ldb = lda; 1445 PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb)); 1446 /* Ainv is Q, overwritten with inverse */ 1447 1448 #if defined(PETSC_USE_COMPLEX) 1449 { 1450 PetscInt i; 1451 for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]); 1452 ierr = PetscFree2(A,Ainv);CHKERRQ(ierr); 1453 } 1454 #endif 1455 PetscFunctionReturn(0); 1456 } 1457 1458 /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */ 1459 static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B) 1460 { 1461 PetscErrorCode ierr; 1462 PetscReal *Bv; 1463 PetscInt i,j; 1464 1465 PetscFunctionBegin; 1466 ierr = PetscMalloc1((ninterval+1)*ndegree,&Bv);CHKERRQ(ierr); 1467 /* Point evaluation of L_p on all the source vertices */ 1468 ierr = PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);CHKERRQ(ierr); 1469 /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */ 1470 for (i=0; i<ninterval; i++) { 1471 for (j=0; j<ndegree; j++) { 1472 if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j]; 1473 else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j]; 1474 } 1475 } 1476 ierr = PetscFree(Bv);CHKERRQ(ierr); 1477 PetscFunctionReturn(0); 1478 } 1479 1480 /*@ 1481 PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals 1482 1483 Not Collective 1484 1485 Input Arguments: 1486 + degree - degree of reconstruction polynomial 1487 . nsource - number of source intervals 1488 . sourcex - sorted coordinates of source cell boundaries (length nsource+1) 1489 . ntarget - number of target intervals 1490 - targetx - sorted coordinates of target cell boundaries (length ntarget+1) 1491 1492 Output Arguments: 1493 . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s] 1494 1495 Level: advanced 1496 1497 .seealso: PetscDTLegendreEval() 1498 @*/ 1499 PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R) 1500 { 1501 PetscErrorCode ierr; 1502 PetscInt i,j,k,*bdegrees,worksize; 1503 PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget; 1504 PetscScalar *tau,*work; 1505 1506 PetscFunctionBegin; 1507 PetscValidRealPointer(sourcex,3); 1508 PetscValidRealPointer(targetx,5); 1509 PetscValidRealPointer(R,6); 1510 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); 1511 #if defined(PETSC_USE_DEBUG) 1512 for (i=0; i<nsource; i++) { 1513 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]); 1514 } 1515 for (i=0; i<ntarget; i++) { 1516 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]); 1517 } 1518 #endif 1519 xmin = PetscMin(sourcex[0],targetx[0]); 1520 xmax = PetscMax(sourcex[nsource],targetx[ntarget]); 1521 center = (xmin + xmax)/2; 1522 hscale = (xmax - xmin)/2; 1523 worksize = nsource; 1524 ierr = PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);CHKERRQ(ierr); 1525 ierr = PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);CHKERRQ(ierr); 1526 for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale; 1527 for (i=0; i<=degree; i++) bdegrees[i] = i+1; 1528 ierr = PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);CHKERRQ(ierr); 1529 ierr = PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);CHKERRQ(ierr); 1530 for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale; 1531 ierr = PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);CHKERRQ(ierr); 1532 for (i=0; i<ntarget; i++) { 1533 PetscReal rowsum = 0; 1534 for (j=0; j<nsource; j++) { 1535 PetscReal sum = 0; 1536 for (k=0; k<degree+1; k++) { 1537 sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j]; 1538 } 1539 R[i*nsource+j] = sum; 1540 rowsum += sum; 1541 } 1542 for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */ 1543 } 1544 ierr = PetscFree4(bdegrees,sourcey,Bsource,work);CHKERRQ(ierr); 1545 ierr = PetscFree4(tau,Bsinv,targety,Btarget);CHKERRQ(ierr); 1546 PetscFunctionReturn(0); 1547 } 1548 1549 /*@C 1550 PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points 1551 1552 Not Collective 1553 1554 Input Parameter: 1555 + n - the number of GLL nodes 1556 . nodes - the GLL nodes 1557 . weights - the GLL weights 1558 . f - the function values at the nodes 1559 1560 Output Parameter: 1561 . in - the value of the integral 1562 1563 Level: beginner 1564 1565 .seealso: PetscDTGaussLobattoLegendreQuadrature() 1566 1567 @*/ 1568 PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal *nodes,PetscReal *weights,const PetscReal *f,PetscReal *in) 1569 { 1570 PetscInt i; 1571 1572 PetscFunctionBegin; 1573 *in = 0.; 1574 for (i=0; i<n; i++) { 1575 *in += f[i]*f[i]*weights[i]; 1576 } 1577 PetscFunctionReturn(0); 1578 } 1579 1580 /*@C 1581 PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element 1582 1583 Not Collective 1584 1585 Input Parameter: 1586 + n - the number of GLL nodes 1587 . nodes - the GLL nodes 1588 . weights - the GLL weights 1589 1590 Output Parameter: 1591 . A - the stiffness element 1592 1593 Level: beginner 1594 1595 Notes: 1596 Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy() 1597 1598 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) 1599 1600 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy() 1601 1602 @*/ 1603 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1604 { 1605 PetscReal **A; 1606 PetscErrorCode ierr; 1607 const PetscReal *gllnodes = nodes; 1608 const PetscInt p = n-1; 1609 PetscReal z0,z1,z2 = -1,x,Lpj,Lpr; 1610 PetscInt i,j,nn,r; 1611 1612 PetscFunctionBegin; 1613 ierr = PetscMalloc1(n,&A);CHKERRQ(ierr); 1614 ierr = PetscMalloc1(n*n,&A[0]);CHKERRQ(ierr); 1615 for (i=1; i<n; i++) A[i] = A[i-1]+n; 1616 1617 for (j=1; j<p; j++) { 1618 x = gllnodes[j]; 1619 z0 = 1.; 1620 z1 = x; 1621 for (nn=1; nn<p; nn++) { 1622 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1623 z0 = z1; 1624 z1 = z2; 1625 } 1626 Lpj=z2; 1627 for (r=1; r<p; r++) { 1628 if (r == j) { 1629 A[j][j]=2./(3.*(1.-gllnodes[j]*gllnodes[j])*Lpj*Lpj); 1630 } else { 1631 x = gllnodes[r]; 1632 z0 = 1.; 1633 z1 = x; 1634 for (nn=1; nn<p; nn++) { 1635 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1636 z0 = z1; 1637 z1 = z2; 1638 } 1639 Lpr = z2; 1640 A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(gllnodes[j]-gllnodes[r])*(gllnodes[j]-gllnodes[r])); 1641 } 1642 } 1643 } 1644 for (j=1; j<p+1; j++) { 1645 x = gllnodes[j]; 1646 z0 = 1.; 1647 z1 = x; 1648 for (nn=1; nn<p; nn++) { 1649 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1650 z0 = z1; 1651 z1 = z2; 1652 } 1653 Lpj = z2; 1654 A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+gllnodes[j])*(1.+gllnodes[j])); 1655 A[0][j] = A[j][0]; 1656 } 1657 for (j=0; j<p; j++) { 1658 x = gllnodes[j]; 1659 z0 = 1.; 1660 z1 = x; 1661 for (nn=1; nn<p; nn++) { 1662 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1663 z0 = z1; 1664 z1 = z2; 1665 } 1666 Lpj=z2; 1667 1668 A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-gllnodes[j])*(1.-gllnodes[j])); 1669 A[j][p] = A[p][j]; 1670 } 1671 A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.; 1672 A[p][p]=A[0][0]; 1673 *AA = A; 1674 PetscFunctionReturn(0); 1675 } 1676 1677 /*@C 1678 PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element 1679 1680 Not Collective 1681 1682 Input Parameter: 1683 + n - the number of GLL nodes 1684 . nodes - the GLL nodes 1685 . weights - the GLL weightss 1686 - A - the stiffness element 1687 1688 Level: beginner 1689 1690 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate() 1691 1692 @*/ 1693 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1694 { 1695 PetscErrorCode ierr; 1696 1697 PetscFunctionBegin; 1698 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 1699 ierr = PetscFree(*AA);CHKERRQ(ierr); 1700 *AA = NULL; 1701 PetscFunctionReturn(0); 1702 } 1703 1704 /*@C 1705 PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element 1706 1707 Not Collective 1708 1709 Input Parameter: 1710 + n - the number of GLL nodes 1711 . nodes - the GLL nodes 1712 . weights - the GLL weights 1713 1714 Output Parameter: 1715 . AA - the stiffness element 1716 - AAT - the transpose of AA (pass in NULL if you do not need this array) 1717 1718 Level: beginner 1719 1720 Notes: 1721 Destroy this with PetscGaussLobattoLegendreElementGradientDestroy() 1722 1723 You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented 1724 1725 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy() 1726 1727 @*/ 1728 PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT) 1729 { 1730 PetscReal **A, **AT = NULL; 1731 PetscErrorCode ierr; 1732 const PetscReal *gllnodes = nodes; 1733 const PetscInt p = n-1; 1734 PetscReal q,qp,Li, Lj,d0; 1735 PetscInt i,j; 1736 1737 PetscFunctionBegin; 1738 ierr = PetscMalloc1(n,&A);CHKERRQ(ierr); 1739 ierr = PetscMalloc1(n*n,&A[0]);CHKERRQ(ierr); 1740 for (i=1; i<n; i++) A[i] = A[i-1]+n; 1741 1742 if (AAT) { 1743 ierr = PetscMalloc1(n,&AT);CHKERRQ(ierr); 1744 ierr = PetscMalloc1(n*n,&AT[0]);CHKERRQ(ierr); 1745 for (i=1; i<n; i++) AT[i] = AT[i-1]+n; 1746 } 1747 1748 if (n==1) {A[0][0] = 0.;} 1749 d0 = (PetscReal)p*((PetscReal)p+1.)/4.; 1750 for (i=0; i<n; i++) { 1751 for (j=0; j<n; j++) { 1752 A[i][j] = 0.; 1753 qAndLEvaluation(p,gllnodes[i],&q,&qp,&Li); 1754 qAndLEvaluation(p,gllnodes[j],&q,&qp,&Lj); 1755 if (i!=j) A[i][j] = Li/(Lj*(gllnodes[i]-gllnodes[j])); 1756 if ((j==i) && (i==0)) A[i][j] = -d0; 1757 if (j==i && i==p) A[i][j] = d0; 1758 if (AT) AT[j][i] = A[i][j]; 1759 } 1760 } 1761 if (AAT) *AAT = AT; 1762 *AA = A; 1763 PetscFunctionReturn(0); 1764 } 1765 1766 /*@C 1767 PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGaussLobattoLegendreElementGradientCreate() 1768 1769 Not Collective 1770 1771 Input Parameter: 1772 + n - the number of GLL nodes 1773 . nodes - the GLL nodes 1774 . weights - the GLL weights 1775 . AA - the stiffness element 1776 - AAT - the transpose of the element 1777 1778 Level: beginner 1779 1780 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionCreate() 1781 1782 @*/ 1783 PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT) 1784 { 1785 PetscErrorCode ierr; 1786 1787 PetscFunctionBegin; 1788 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 1789 ierr = PetscFree(*AA);CHKERRQ(ierr); 1790 *AA = NULL; 1791 if (*AAT) { 1792 ierr = PetscFree((*AAT)[0]);CHKERRQ(ierr); 1793 ierr = PetscFree(*AAT);CHKERRQ(ierr); 1794 *AAT = NULL; 1795 } 1796 PetscFunctionReturn(0); 1797 } 1798 1799 /*@C 1800 PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element 1801 1802 Not Collective 1803 1804 Input Parameter: 1805 + n - the number of GLL nodes 1806 . nodes - the GLL nodes 1807 . weights - the GLL weightss 1808 1809 Output Parameter: 1810 . AA - the stiffness element 1811 1812 Level: beginner 1813 1814 Notes: 1815 Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy() 1816 1817 This is the same as the Gradient operator multiplied by the diagonal mass matrix 1818 1819 You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented 1820 1821 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionDestroy() 1822 1823 @*/ 1824 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1825 { 1826 PetscReal **D; 1827 PetscErrorCode ierr; 1828 const PetscReal *gllweights = weights; 1829 const PetscInt glln = n; 1830 PetscInt i,j; 1831 1832 PetscFunctionBegin; 1833 ierr = PetscGaussLobattoLegendreElementGradientCreate(n,nodes,weights,&D,NULL);CHKERRQ(ierr); 1834 for (i=0; i<glln; i++){ 1835 for (j=0; j<glln; j++) { 1836 D[i][j] = gllweights[i]*D[i][j]; 1837 } 1838 } 1839 *AA = D; 1840 PetscFunctionReturn(0); 1841 } 1842 1843 /*@C 1844 PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element 1845 1846 Not Collective 1847 1848 Input Parameter: 1849 + n - the number of GLL nodes 1850 . nodes - the GLL nodes 1851 . weights - the GLL weights 1852 - A - advection 1853 1854 Level: beginner 1855 1856 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementAdvectionCreate() 1857 1858 @*/ 1859 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1860 { 1861 PetscErrorCode ierr; 1862 1863 PetscFunctionBegin; 1864 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 1865 ierr = PetscFree(*AA);CHKERRQ(ierr); 1866 *AA = NULL; 1867 PetscFunctionReturn(0); 1868 } 1869 1870 PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1871 { 1872 PetscReal **A; 1873 PetscErrorCode ierr; 1874 const PetscReal *gllweights = weights; 1875 const PetscInt glln = n; 1876 PetscInt i,j; 1877 1878 PetscFunctionBegin; 1879 ierr = PetscMalloc1(glln,&A);CHKERRQ(ierr); 1880 ierr = PetscMalloc1(glln*glln,&A[0]);CHKERRQ(ierr); 1881 for (i=1; i<glln; i++) A[i] = A[i-1]+glln; 1882 if (glln==1) {A[0][0] = 0.;} 1883 for (i=0; i<glln; i++) { 1884 for (j=0; j<glln; j++) { 1885 A[i][j] = 0.; 1886 if (j==i) A[i][j] = gllweights[i]; 1887 } 1888 } 1889 *AA = A; 1890 PetscFunctionReturn(0); 1891 } 1892 1893 PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1894 { 1895 PetscErrorCode ierr; 1896 1897 PetscFunctionBegin; 1898 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 1899 ierr = PetscFree(*AA);CHKERRQ(ierr); 1900 *AA = NULL; 1901 PetscFunctionReturn(0); 1902 } 1903