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 PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info)); 697 ierr = PetscFPTrapPop();CHKERRQ(ierr); 698 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error"); 699 700 for (i=0; i<(npoints+1)/2; i++) { 701 PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */ 702 x[i] = (a+b)/2 - y*(b-a)/2; 703 if (x[i] == -0.0) x[i] = 0.0; 704 x[npoints-i-1] = (a+b)/2 + y*(b-a)/2; 705 706 w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints]))); 707 } 708 ierr = PetscFree2(Z,work);CHKERRQ(ierr); 709 PetscFunctionReturn(0); 710 } 711 712 static void qAndLEvaluation(PetscInt n, PetscReal x, PetscReal *q, PetscReal *qp, PetscReal *Ln) 713 /* 714 Compute the polynomial q(x) = L_{N+1}(x) - L_{n-1}(x) and its derivative in 715 addition to L_N(x) as these are needed for computing the GLL points via Newton's method. 716 Reference: "Implementing Spectral Methods for Partial Differential Equations: Algorithms 717 for Scientists and Engineers" by David A. Kopriva. 718 */ 719 { 720 PetscInt k; 721 722 PetscReal Lnp; 723 PetscReal Lnp1, Lnp1p; 724 PetscReal Lnm1, Lnm1p; 725 PetscReal Lnm2, Lnm2p; 726 727 Lnm1 = 1.0; 728 *Ln = x; 729 Lnm1p = 0.0; 730 Lnp = 1.0; 731 732 for (k=2; k<=n; ++k) { 733 Lnm2 = Lnm1; 734 Lnm1 = *Ln; 735 Lnm2p = Lnm1p; 736 Lnm1p = Lnp; 737 *Ln = (2.*((PetscReal)k)-1.)/(1.0*((PetscReal)k))*x*Lnm1 - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm2; 738 Lnp = Lnm2p + (2.0*((PetscReal)k)-1.)*Lnm1; 739 } 740 k = n+1; 741 Lnp1 = (2.*((PetscReal)k)-1.)/(((PetscReal)k))*x*(*Ln) - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm1; 742 Lnp1p = Lnm1p + (2.0*((PetscReal)k)-1.)*(*Ln); 743 *q = Lnp1 - Lnm1; 744 *qp = Lnp1p - Lnm1p; 745 } 746 747 /*@C 748 PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre 749 nodes of a given size on the domain [-1,1] 750 751 Not Collective 752 753 Input Parameter: 754 + n - number of grid nodes 755 - type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON 756 757 Output Arguments: 758 + x - quadrature points 759 - w - quadrature weights 760 761 Notes: 762 For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not 763 close enough to the desired solution 764 765 These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes 766 767 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 768 769 Level: intermediate 770 771 .seealso: PetscDTGaussQuadrature() 772 773 @*/ 774 PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal *x,PetscReal *w) 775 { 776 PetscErrorCode ierr; 777 778 PetscFunctionBegin; 779 if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide at least 2 grid points per element"); 780 781 if (type == PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA) { 782 PetscReal *M,si; 783 PetscBLASInt bn,lierr; 784 PetscReal x0,z0,z1,z2; 785 PetscInt i,p = npoints - 1,nn; 786 787 x[0] =-1.0; 788 x[npoints-1] = 1.0; 789 if (npoints-2 > 0){ 790 ierr = PetscMalloc1(npoints-1,&M);CHKERRQ(ierr); 791 for (i=0; i<npoints-2; i++) { 792 si = ((PetscReal)i)+1.0; 793 M[i]=0.5*PetscSqrtReal(si*(si+2.0)/((si+0.5)*(si+1.5))); 794 } 795 ierr = PetscBLASIntCast(npoints-2,&bn);CHKERRQ(ierr); 796 ierr = PetscArrayzero(&x[1],bn);CHKERRQ(ierr); 797 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 798 x0=0; 799 PetscStackCallBLAS("LAPACKsteqr",LAPACKREALsteqr_("N",&bn,&x[1],M,&x0,&bn,M,&lierr)); 800 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in STERF Lapack routine %d",(int)lierr); 801 ierr = PetscFPTrapPop();CHKERRQ(ierr); 802 ierr = PetscFree(M);CHKERRQ(ierr); 803 } 804 if ((npoints-1)%2==0) { 805 x[(npoints-1)/2] = 0.0; /* hard wire to exactly 0.0 since linear algebra produces nonzero */ 806 } 807 808 w[0] = w[p] = 2.0/(((PetscReal)(p))*(((PetscReal)p)+1.0)); 809 z2 = -1.; /* Dummy value to avoid -Wmaybe-initialized */ 810 for (i=1; i<p; i++) { 811 x0 = x[i]; 812 z0 = 1.0; 813 z1 = x0; 814 for (nn=1; nn<p; nn++) { 815 z2 = x0*z1*(2.0*((PetscReal)nn)+1.0)/(((PetscReal)nn)+1.0)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.0)); 816 z0 = z1; 817 z1 = z2; 818 } 819 w[i]=2.0/(((PetscReal)p)*(((PetscReal)p)+1.0)*z2*z2); 820 } 821 } else { 822 PetscInt j,m; 823 PetscReal z1,z,q,qp,Ln; 824 PetscReal *pt; 825 ierr = PetscMalloc1(npoints,&pt);CHKERRQ(ierr); 826 827 if (npoints > 30) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON produces incorrect answers for n > 30"); 828 x[0] = -1.0; 829 x[npoints-1] = 1.0; 830 w[0] = w[npoints-1] = 2./(((PetscReal)npoints)*(((PetscReal)npoints)-1.0)); 831 m = (npoints-1)/2; /* The roots are symmetric, so we only find half of them. */ 832 for (j=1; j<=m; j++) { /* Loop over the desired roots. */ 833 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)); 834 /* Starting with the above approximation to the ith root, we enter */ 835 /* the main loop of refinement by Newton's method. */ 836 do { 837 qAndLEvaluation(npoints-1,z,&q,&qp,&Ln); 838 z1 = z; 839 z = z1-q/qp; /* Newton's method. */ 840 } while (PetscAbs(z-z1) > 10.*PETSC_MACHINE_EPSILON); 841 qAndLEvaluation(npoints-1,z,&q,&qp,&Ln); 842 843 x[j] = z; 844 x[npoints-1-j] = -z; /* and put in its symmetric counterpart. */ 845 w[j] = 2.0/(((PetscReal)npoints)*(((PetscReal)npoints)-1.)*Ln*Ln); /* Compute the weight */ 846 w[npoints-1-j] = w[j]; /* and its symmetric counterpart. */ 847 pt[j]=qp; 848 } 849 850 if ((npoints-1)%2==0) { 851 qAndLEvaluation(npoints-1,0.0,&q,&qp,&Ln); 852 x[(npoints-1)/2] = 0.0; 853 w[(npoints-1)/2] = 2.0/(((PetscReal)npoints)*(((PetscReal)npoints)-1.)*Ln*Ln); 854 } 855 ierr = PetscFree(pt);CHKERRQ(ierr); 856 } 857 PetscFunctionReturn(0); 858 } 859 860 /*@ 861 PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature 862 863 Not Collective 864 865 Input Arguments: 866 + dim - The spatial dimension 867 . Nc - The number of components 868 . npoints - number of points in one dimension 869 . a - left end of interval (often-1) 870 - b - right end of interval (often +1) 871 872 Output Argument: 873 . q - A PetscQuadrature object 874 875 Level: intermediate 876 877 .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval() 878 @*/ 879 PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 880 { 881 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c; 882 PetscReal *x, *w, *xw, *ww; 883 PetscErrorCode ierr; 884 885 PetscFunctionBegin; 886 ierr = PetscMalloc1(totpoints*dim,&x);CHKERRQ(ierr); 887 ierr = PetscMalloc1(totpoints*Nc,&w);CHKERRQ(ierr); 888 /* Set up the Golub-Welsch system */ 889 switch (dim) { 890 case 0: 891 ierr = PetscFree(x);CHKERRQ(ierr); 892 ierr = PetscFree(w);CHKERRQ(ierr); 893 ierr = PetscMalloc1(1, &x);CHKERRQ(ierr); 894 ierr = PetscMalloc1(Nc, &w);CHKERRQ(ierr); 895 x[0] = 0.0; 896 for (c = 0; c < Nc; ++c) w[c] = 1.0; 897 break; 898 case 1: 899 ierr = PetscMalloc1(npoints,&ww);CHKERRQ(ierr); 900 ierr = PetscDTGaussQuadrature(npoints, a, b, x, ww);CHKERRQ(ierr); 901 for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i]; 902 ierr = PetscFree(ww);CHKERRQ(ierr); 903 break; 904 case 2: 905 ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr); 906 ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr); 907 for (i = 0; i < npoints; ++i) { 908 for (j = 0; j < npoints; ++j) { 909 x[(i*npoints+j)*dim+0] = xw[i]; 910 x[(i*npoints+j)*dim+1] = xw[j]; 911 for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j]; 912 } 913 } 914 ierr = PetscFree2(xw,ww);CHKERRQ(ierr); 915 break; 916 case 3: 917 ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr); 918 ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr); 919 for (i = 0; i < npoints; ++i) { 920 for (j = 0; j < npoints; ++j) { 921 for (k = 0; k < npoints; ++k) { 922 x[((i*npoints+j)*npoints+k)*dim+0] = xw[i]; 923 x[((i*npoints+j)*npoints+k)*dim+1] = xw[j]; 924 x[((i*npoints+j)*npoints+k)*dim+2] = xw[k]; 925 for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k]; 926 } 927 } 928 } 929 ierr = PetscFree2(xw,ww);CHKERRQ(ierr); 930 break; 931 default: 932 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim); 933 } 934 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 935 ierr = PetscQuadratureSetOrder(*q, 2*npoints-1);CHKERRQ(ierr); 936 ierr = PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);CHKERRQ(ierr); 937 ierr = PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");CHKERRQ(ierr); 938 PetscFunctionReturn(0); 939 } 940 941 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x. 942 Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */ 943 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 944 { 945 PetscReal apb, pn1, pn2; 946 PetscInt k; 947 948 PetscFunctionBegin; 949 if (!n) {*P = 1.0; PetscFunctionReturn(0);} 950 if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); PetscFunctionReturn(0);} 951 apb = a + b; 952 pn2 = 1.0; 953 pn1 = 0.5 * (a - b + (apb + 2.0) * x); 954 *P = 0.0; 955 for (k = 2; k < n+1; ++k) { 956 PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0); 957 PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b); 958 PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb); 959 PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb); 960 961 a2 = a2 / a1; 962 a3 = a3 / a1; 963 a4 = a4 / a1; 964 *P = (a2 + a3 * x) * pn1 - a4 * pn2; 965 pn2 = pn1; 966 pn1 = *P; 967 } 968 PetscFunctionReturn(0); 969 } 970 971 /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */ 972 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 973 { 974 PetscReal nP; 975 PetscErrorCode ierr; 976 977 PetscFunctionBegin; 978 if (!n) {*P = 0.0; PetscFunctionReturn(0);} 979 ierr = PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);CHKERRQ(ierr); 980 *P = 0.5 * (a + b + n + 1) * nP; 981 PetscFunctionReturn(0); 982 } 983 984 /* Maps from [-1,1]^2 to the (-1,1) reference triangle */ 985 PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta) 986 { 987 PetscFunctionBegin; 988 *xi = 0.5 * (1.0 + x) * (1.0 - y) - 1.0; 989 *eta = y; 990 PetscFunctionReturn(0); 991 } 992 993 /* Maps from [-1,1]^2 to the (-1,1) reference triangle */ 994 PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta) 995 { 996 PetscFunctionBegin; 997 *xi = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0; 998 *eta = 0.5 * (1.0 + y) * (1.0 - z) - 1.0; 999 *zeta = z; 1000 PetscFunctionReturn(0); 1001 } 1002 1003 static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w) 1004 { 1005 PetscInt maxIter = 100; 1006 PetscReal eps = 1.0e-8; 1007 PetscReal a1, a2, a3, a4, a5, a6; 1008 PetscInt k; 1009 PetscErrorCode ierr; 1010 1011 PetscFunctionBegin; 1012 1013 a1 = PetscPowReal(2.0, a+b+1); 1014 #if defined(PETSC_HAVE_TGAMMA) 1015 a2 = PetscTGamma(a + npoints + 1); 1016 a3 = PetscTGamma(b + npoints + 1); 1017 a4 = PetscTGamma(a + b + npoints + 1); 1018 #else 1019 { 1020 PetscInt ia, ib; 1021 1022 ia = (PetscInt) a; 1023 ib = (PetscInt) b; 1024 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 */ 1025 ierr = PetscDTFactorial(ia + npoints, &a2);CHKERRQ(ierr); 1026 ierr = PetscDTFactorial(ib + npoints, &a3);CHKERRQ(ierr); 1027 ierr = PetscDTFactorial(ia + ib + npoints, &a4);CHKERRQ(ierr); 1028 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable."); 1029 } 1030 #endif 1031 1032 ierr = PetscDTFactorial(npoints, &a5);CHKERRQ(ierr); 1033 a6 = a1 * a2 * a3 / a4 / a5; 1034 /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses. 1035 Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */ 1036 for (k = 0; k < npoints; ++k) { 1037 PetscReal r = -PetscCosReal((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP; 1038 PetscInt j; 1039 1040 if (k > 0) r = 0.5 * (r + x[k-1]); 1041 for (j = 0; j < maxIter; ++j) { 1042 PetscReal s = 0.0, delta, f, fp; 1043 PetscInt i; 1044 1045 for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]); 1046 ierr = PetscDTComputeJacobi(a, b, npoints, r, &f);CHKERRQ(ierr); 1047 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);CHKERRQ(ierr); 1048 delta = f / (fp - f * s); 1049 r = r - delta; 1050 if (PetscAbsReal(delta) < eps) break; 1051 } 1052 x[k] = r; 1053 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);CHKERRQ(ierr); 1054 w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP); 1055 } 1056 PetscFunctionReturn(0); 1057 } 1058 1059 /*@ 1060 PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex 1061 1062 Not Collective 1063 1064 Input Arguments: 1065 + dim - The simplex dimension 1066 . Nc - The number of components 1067 . npoints - The number of points in one dimension 1068 . a - left end of interval (often-1) 1069 - b - right end of interval (often +1) 1070 1071 Output Argument: 1072 . q - A PetscQuadrature object 1073 1074 Level: intermediate 1075 1076 References: 1077 . 1. - Karniadakis and Sherwin. FIAT 1078 1079 .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature() 1080 @*/ 1081 PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q) 1082 { 1083 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints; 1084 PetscReal *px, *wx, *py, *wy, *pz, *wz, *x, *w; 1085 PetscInt i, j, k, c; 1086 PetscErrorCode ierr; 1087 1088 PetscFunctionBegin; 1089 if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now"); 1090 ierr = PetscMalloc1(totpoints*dim, &x);CHKERRQ(ierr); 1091 ierr = PetscMalloc1(totpoints*Nc, &w);CHKERRQ(ierr); 1092 switch (dim) { 1093 case 0: 1094 ierr = PetscFree(x);CHKERRQ(ierr); 1095 ierr = PetscFree(w);CHKERRQ(ierr); 1096 ierr = PetscMalloc1(1, &x);CHKERRQ(ierr); 1097 ierr = PetscMalloc1(Nc, &w);CHKERRQ(ierr); 1098 x[0] = 0.0; 1099 for (c = 0; c < Nc; ++c) w[c] = 1.0; 1100 break; 1101 case 1: 1102 ierr = PetscMalloc1(npoints,&wx);CHKERRQ(ierr); 1103 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, wx);CHKERRQ(ierr); 1104 for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i]; 1105 ierr = PetscFree(wx);CHKERRQ(ierr); 1106 break; 1107 case 2: 1108 ierr = PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);CHKERRQ(ierr); 1109 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);CHKERRQ(ierr); 1110 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);CHKERRQ(ierr); 1111 for (i = 0; i < npoints; ++i) { 1112 for (j = 0; j < npoints; ++j) { 1113 ierr = PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);CHKERRQ(ierr); 1114 for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j]; 1115 } 1116 } 1117 ierr = PetscFree4(px,wx,py,wy);CHKERRQ(ierr); 1118 break; 1119 case 3: 1120 ierr = PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);CHKERRQ(ierr); 1121 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);CHKERRQ(ierr); 1122 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);CHKERRQ(ierr); 1123 ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);CHKERRQ(ierr); 1124 for (i = 0; i < npoints; ++i) { 1125 for (j = 0; j < npoints; ++j) { 1126 for (k = 0; k < npoints; ++k) { 1127 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); 1128 for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k]; 1129 } 1130 } 1131 } 1132 ierr = PetscFree6(px,wx,py,wy,pz,wz);CHKERRQ(ierr); 1133 break; 1134 default: 1135 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim); 1136 } 1137 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 1138 ierr = PetscQuadratureSetOrder(*q, 2*npoints-1);CHKERRQ(ierr); 1139 ierr = PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);CHKERRQ(ierr); 1140 ierr = PetscObjectChangeTypeName((PetscObject)*q,"GaussJacobi");CHKERRQ(ierr); 1141 PetscFunctionReturn(0); 1142 } 1143 1144 /*@ 1145 PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell 1146 1147 Not Collective 1148 1149 Input Arguments: 1150 + dim - The cell dimension 1151 . level - The number of points in one dimension, 2^l 1152 . a - left end of interval (often-1) 1153 - b - right end of interval (often +1) 1154 1155 Output Argument: 1156 . q - A PetscQuadrature object 1157 1158 Level: intermediate 1159 1160 .seealso: PetscDTGaussTensorQuadrature() 1161 @*/ 1162 PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q) 1163 { 1164 const PetscInt p = 16; /* Digits of precision in the evaluation */ 1165 const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */ 1166 const PetscReal beta = (b+a)/2.; /* Center of the integration interval */ 1167 const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */ 1168 PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */ 1169 PetscReal wk = 0.5*PETSC_PI; /* Quadrature weight at x_k */ 1170 PetscReal *x, *w; 1171 PetscInt K, k, npoints; 1172 PetscErrorCode ierr; 1173 1174 PetscFunctionBegin; 1175 if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim); 1176 if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits"); 1177 /* Find K such that the weights are < 32 digits of precision */ 1178 for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) { 1179 wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h))); 1180 } 1181 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr); 1182 ierr = PetscQuadratureSetOrder(*q, 2*K+1);CHKERRQ(ierr); 1183 npoints = 2*K-1; 1184 ierr = PetscMalloc1(npoints*dim, &x);CHKERRQ(ierr); 1185 ierr = PetscMalloc1(npoints, &w);CHKERRQ(ierr); 1186 /* Center term */ 1187 x[0] = beta; 1188 w[0] = 0.5*alpha*PETSC_PI; 1189 for (k = 1; k < K; ++k) { 1190 wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1191 xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h)); 1192 x[2*k-1] = -alpha*xk+beta; 1193 w[2*k-1] = wk; 1194 x[2*k+0] = alpha*xk+beta; 1195 w[2*k+0] = wk; 1196 } 1197 ierr = PetscQuadratureSetData(*q, dim, 1, npoints, x, w);CHKERRQ(ierr); 1198 PetscFunctionReturn(0); 1199 } 1200 1201 PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 1202 { 1203 const PetscInt p = 16; /* Digits of precision in the evaluation */ 1204 const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */ 1205 const PetscReal beta = (b+a)/2.; /* Center of the integration interval */ 1206 PetscReal h = 1.0; /* Step size, length between x_k */ 1207 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 1208 PetscReal osum = 0.0; /* Integral on last level */ 1209 PetscReal psum = 0.0; /* Integral on the level before the last level */ 1210 PetscReal sum; /* Integral on current level */ 1211 PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 1212 PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 1213 PetscReal wk; /* Quadrature weight at x_k */ 1214 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 1215 PetscInt d; /* Digits of precision in the integral */ 1216 1217 PetscFunctionBegin; 1218 if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 1219 /* Center term */ 1220 func(beta, &lval); 1221 sum = 0.5*alpha*PETSC_PI*lval; 1222 /* */ 1223 do { 1224 PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4; 1225 PetscInt k = 1; 1226 1227 ++l; 1228 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */ 1229 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 1230 psum = osum; 1231 osum = sum; 1232 h *= 0.5; 1233 sum *= 0.5; 1234 do { 1235 wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1236 yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h))); 1237 lx = -alpha*(1.0 - yk)+beta; 1238 rx = alpha*(1.0 - yk)+beta; 1239 func(lx, &lval); 1240 func(rx, &rval); 1241 lterm = alpha*wk*lval; 1242 maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm); 1243 sum += lterm; 1244 rterm = alpha*wk*rval; 1245 maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm); 1246 sum += rterm; 1247 ++k; 1248 /* Only need to evaluate every other point on refined levels */ 1249 if (l != 1) ++k; 1250 } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */ 1251 1252 d1 = PetscLog10Real(PetscAbsReal(sum - osum)); 1253 d2 = PetscLog10Real(PetscAbsReal(sum - psum)); 1254 d3 = PetscLog10Real(maxTerm) - p; 1255 if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0; 1256 else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm))); 1257 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4))); 1258 } while (d < digits && l < 12); 1259 *sol = sum; 1260 1261 PetscFunctionReturn(0); 1262 } 1263 1264 #if defined(PETSC_HAVE_MPFR) 1265 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 1266 { 1267 const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */ 1268 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */ 1269 mpfr_t alpha; /* Half-width of the integration interval */ 1270 mpfr_t beta; /* Center of the integration interval */ 1271 mpfr_t h; /* Step size, length between x_k */ 1272 mpfr_t osum; /* Integral on last level */ 1273 mpfr_t psum; /* Integral on the level before the last level */ 1274 mpfr_t sum; /* Integral on current level */ 1275 mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */ 1276 mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */ 1277 mpfr_t wk; /* Quadrature weight at x_k */ 1278 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */ 1279 PetscInt d; /* Digits of precision in the integral */ 1280 mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp; 1281 1282 PetscFunctionBegin; 1283 if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits"); 1284 /* Create high precision storage */ 1285 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); 1286 /* Initialization */ 1287 mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN); 1288 mpfr_set_d(beta, 0.5*(b+a), MPFR_RNDN); 1289 mpfr_set_d(osum, 0.0, MPFR_RNDN); 1290 mpfr_set_d(psum, 0.0, MPFR_RNDN); 1291 mpfr_set_d(h, 1.0, MPFR_RNDN); 1292 mpfr_const_pi(pi2, MPFR_RNDN); 1293 mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN); 1294 /* Center term */ 1295 func(0.5*(b+a), &lval); 1296 mpfr_set(sum, pi2, MPFR_RNDN); 1297 mpfr_mul(sum, sum, alpha, MPFR_RNDN); 1298 mpfr_mul_d(sum, sum, lval, MPFR_RNDN); 1299 /* */ 1300 do { 1301 PetscReal d1, d2, d3, d4; 1302 PetscInt k = 1; 1303 1304 ++l; 1305 mpfr_set_d(maxTerm, 0.0, MPFR_RNDN); 1306 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */ 1307 /* At each level of refinement, h --> h/2 and sum --> sum/2 */ 1308 mpfr_set(psum, osum, MPFR_RNDN); 1309 mpfr_set(osum, sum, MPFR_RNDN); 1310 mpfr_mul_d(h, h, 0.5, MPFR_RNDN); 1311 mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN); 1312 do { 1313 mpfr_set_si(kh, k, MPFR_RNDN); 1314 mpfr_mul(kh, kh, h, MPFR_RNDN); 1315 /* Weight */ 1316 mpfr_set(wk, h, MPFR_RNDN); 1317 mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN); 1318 mpfr_mul(msinh, msinh, pi2, MPFR_RNDN); 1319 mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN); 1320 mpfr_cosh(tmp, msinh, MPFR_RNDN); 1321 mpfr_sqr(tmp, tmp, MPFR_RNDN); 1322 mpfr_mul(wk, wk, mcosh, MPFR_RNDN); 1323 mpfr_div(wk, wk, tmp, MPFR_RNDN); 1324 /* Abscissa */ 1325 mpfr_set_d(yk, 1.0, MPFR_RNDZ); 1326 mpfr_cosh(tmp, msinh, MPFR_RNDN); 1327 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 1328 mpfr_exp(tmp, msinh, MPFR_RNDN); 1329 mpfr_div(yk, yk, tmp, MPFR_RNDZ); 1330 /* Quadrature points */ 1331 mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ); 1332 mpfr_mul(lx, lx, alpha, MPFR_RNDU); 1333 mpfr_add(lx, lx, beta, MPFR_RNDU); 1334 mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ); 1335 mpfr_mul(rx, rx, alpha, MPFR_RNDD); 1336 mpfr_add(rx, rx, beta, MPFR_RNDD); 1337 /* Evaluation */ 1338 func(mpfr_get_d(lx, MPFR_RNDU), &lval); 1339 func(mpfr_get_d(rx, MPFR_RNDD), &rval); 1340 /* Update */ 1341 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 1342 mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN); 1343 mpfr_add(sum, sum, tmp, MPFR_RNDN); 1344 mpfr_abs(tmp, tmp, MPFR_RNDN); 1345 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 1346 mpfr_set(curTerm, tmp, MPFR_RNDN); 1347 mpfr_mul(tmp, wk, alpha, MPFR_RNDN); 1348 mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN); 1349 mpfr_add(sum, sum, tmp, MPFR_RNDN); 1350 mpfr_abs(tmp, tmp, MPFR_RNDN); 1351 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN); 1352 mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN); 1353 ++k; 1354 /* Only need to evaluate every other point on refined levels */ 1355 if (l != 1) ++k; 1356 mpfr_log10(tmp, wk, MPFR_RNDN); 1357 mpfr_abs(tmp, tmp, MPFR_RNDN); 1358 } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */ 1359 mpfr_sub(tmp, sum, osum, MPFR_RNDN); 1360 mpfr_abs(tmp, tmp, MPFR_RNDN); 1361 mpfr_log10(tmp, tmp, MPFR_RNDN); 1362 d1 = mpfr_get_d(tmp, MPFR_RNDN); 1363 mpfr_sub(tmp, sum, psum, MPFR_RNDN); 1364 mpfr_abs(tmp, tmp, MPFR_RNDN); 1365 mpfr_log10(tmp, tmp, MPFR_RNDN); 1366 d2 = mpfr_get_d(tmp, MPFR_RNDN); 1367 mpfr_log10(tmp, maxTerm, MPFR_RNDN); 1368 d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits; 1369 mpfr_log10(tmp, curTerm, MPFR_RNDN); 1370 d4 = mpfr_get_d(tmp, MPFR_RNDN); 1371 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4))); 1372 } while (d < digits && l < 8); 1373 *sol = mpfr_get_d(sum, MPFR_RNDN); 1374 /* Cleanup */ 1375 mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL); 1376 PetscFunctionReturn(0); 1377 } 1378 #else 1379 1380 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol) 1381 { 1382 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp"); 1383 } 1384 #endif 1385 1386 /* Overwrites A. Can only handle full-rank problems with m>=n 1387 * A in column-major format 1388 * Ainv in row-major format 1389 * tau has length m 1390 * worksize must be >= max(1,n) 1391 */ 1392 static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work) 1393 { 1394 PetscErrorCode ierr; 1395 PetscBLASInt M,N,K,lda,ldb,ldwork,info; 1396 PetscScalar *A,*Ainv,*R,*Q,Alpha; 1397 1398 PetscFunctionBegin; 1399 #if defined(PETSC_USE_COMPLEX) 1400 { 1401 PetscInt i,j; 1402 ierr = PetscMalloc2(m*n,&A,m*n,&Ainv);CHKERRQ(ierr); 1403 for (j=0; j<n; j++) { 1404 for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j]; 1405 } 1406 mstride = m; 1407 } 1408 #else 1409 A = A_in; 1410 Ainv = Ainv_out; 1411 #endif 1412 1413 ierr = PetscBLASIntCast(m,&M);CHKERRQ(ierr); 1414 ierr = PetscBLASIntCast(n,&N);CHKERRQ(ierr); 1415 ierr = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr); 1416 ierr = PetscBLASIntCast(worksize,&ldwork);CHKERRQ(ierr); 1417 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1418 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info)); 1419 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1420 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error"); 1421 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */ 1422 1423 /* Extract an explicit representation of Q */ 1424 Q = Ainv; 1425 ierr = PetscArraycpy(Q,A,mstride*n);CHKERRQ(ierr); 1426 K = N; /* full rank */ 1427 PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info)); 1428 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error"); 1429 1430 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */ 1431 Alpha = 1.0; 1432 ldb = lda; 1433 PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb)); 1434 /* Ainv is Q, overwritten with inverse */ 1435 1436 #if defined(PETSC_USE_COMPLEX) 1437 { 1438 PetscInt i; 1439 for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]); 1440 ierr = PetscFree2(A,Ainv);CHKERRQ(ierr); 1441 } 1442 #endif 1443 PetscFunctionReturn(0); 1444 } 1445 1446 /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */ 1447 static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B) 1448 { 1449 PetscErrorCode ierr; 1450 PetscReal *Bv; 1451 PetscInt i,j; 1452 1453 PetscFunctionBegin; 1454 ierr = PetscMalloc1((ninterval+1)*ndegree,&Bv);CHKERRQ(ierr); 1455 /* Point evaluation of L_p on all the source vertices */ 1456 ierr = PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);CHKERRQ(ierr); 1457 /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */ 1458 for (i=0; i<ninterval; i++) { 1459 for (j=0; j<ndegree; j++) { 1460 if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j]; 1461 else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j]; 1462 } 1463 } 1464 ierr = PetscFree(Bv);CHKERRQ(ierr); 1465 PetscFunctionReturn(0); 1466 } 1467 1468 /*@ 1469 PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals 1470 1471 Not Collective 1472 1473 Input Arguments: 1474 + degree - degree of reconstruction polynomial 1475 . nsource - number of source intervals 1476 . sourcex - sorted coordinates of source cell boundaries (length nsource+1) 1477 . ntarget - number of target intervals 1478 - targetx - sorted coordinates of target cell boundaries (length ntarget+1) 1479 1480 Output Arguments: 1481 . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s] 1482 1483 Level: advanced 1484 1485 .seealso: PetscDTLegendreEval() 1486 @*/ 1487 PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R) 1488 { 1489 PetscErrorCode ierr; 1490 PetscInt i,j,k,*bdegrees,worksize; 1491 PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget; 1492 PetscScalar *tau,*work; 1493 1494 PetscFunctionBegin; 1495 PetscValidRealPointer(sourcex,3); 1496 PetscValidRealPointer(targetx,5); 1497 PetscValidRealPointer(R,6); 1498 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); 1499 #if defined(PETSC_USE_DEBUG) 1500 for (i=0; i<nsource; i++) { 1501 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]); 1502 } 1503 for (i=0; i<ntarget; i++) { 1504 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]); 1505 } 1506 #endif 1507 xmin = PetscMin(sourcex[0],targetx[0]); 1508 xmax = PetscMax(sourcex[nsource],targetx[ntarget]); 1509 center = (xmin + xmax)/2; 1510 hscale = (xmax - xmin)/2; 1511 worksize = nsource; 1512 ierr = PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);CHKERRQ(ierr); 1513 ierr = PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);CHKERRQ(ierr); 1514 for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale; 1515 for (i=0; i<=degree; i++) bdegrees[i] = i+1; 1516 ierr = PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);CHKERRQ(ierr); 1517 ierr = PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);CHKERRQ(ierr); 1518 for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale; 1519 ierr = PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);CHKERRQ(ierr); 1520 for (i=0; i<ntarget; i++) { 1521 PetscReal rowsum = 0; 1522 for (j=0; j<nsource; j++) { 1523 PetscReal sum = 0; 1524 for (k=0; k<degree+1; k++) { 1525 sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j]; 1526 } 1527 R[i*nsource+j] = sum; 1528 rowsum += sum; 1529 } 1530 for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */ 1531 } 1532 ierr = PetscFree4(bdegrees,sourcey,Bsource,work);CHKERRQ(ierr); 1533 ierr = PetscFree4(tau,Bsinv,targety,Btarget);CHKERRQ(ierr); 1534 PetscFunctionReturn(0); 1535 } 1536 1537 /*@C 1538 PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points 1539 1540 Not Collective 1541 1542 Input Parameter: 1543 + n - the number of GLL nodes 1544 . nodes - the GLL nodes 1545 . weights - the GLL weights 1546 . f - the function values at the nodes 1547 1548 Output Parameter: 1549 . in - the value of the integral 1550 1551 Level: beginner 1552 1553 .seealso: PetscDTGaussLobattoLegendreQuadrature() 1554 1555 @*/ 1556 PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal *nodes,PetscReal *weights,const PetscReal *f,PetscReal *in) 1557 { 1558 PetscInt i; 1559 1560 PetscFunctionBegin; 1561 *in = 0.; 1562 for (i=0; i<n; i++) { 1563 *in += f[i]*f[i]*weights[i]; 1564 } 1565 PetscFunctionReturn(0); 1566 } 1567 1568 /*@C 1569 PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element 1570 1571 Not Collective 1572 1573 Input Parameter: 1574 + n - the number of GLL nodes 1575 . nodes - the GLL nodes 1576 . weights - the GLL weights 1577 1578 Output Parameter: 1579 . A - the stiffness element 1580 1581 Level: beginner 1582 1583 Notes: 1584 Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy() 1585 1586 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) 1587 1588 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy() 1589 1590 @*/ 1591 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1592 { 1593 PetscReal **A; 1594 PetscErrorCode ierr; 1595 const PetscReal *gllnodes = nodes; 1596 const PetscInt p = n-1; 1597 PetscReal z0,z1,z2 = -1,x,Lpj,Lpr; 1598 PetscInt i,j,nn,r; 1599 1600 PetscFunctionBegin; 1601 ierr = PetscMalloc1(n,&A);CHKERRQ(ierr); 1602 ierr = PetscMalloc1(n*n,&A[0]);CHKERRQ(ierr); 1603 for (i=1; i<n; i++) A[i] = A[i-1]+n; 1604 1605 for (j=1; j<p; j++) { 1606 x = gllnodes[j]; 1607 z0 = 1.; 1608 z1 = x; 1609 for (nn=1; nn<p; nn++) { 1610 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1611 z0 = z1; 1612 z1 = z2; 1613 } 1614 Lpj=z2; 1615 for (r=1; r<p; r++) { 1616 if (r == j) { 1617 A[j][j]=2./(3.*(1.-gllnodes[j]*gllnodes[j])*Lpj*Lpj); 1618 } else { 1619 x = gllnodes[r]; 1620 z0 = 1.; 1621 z1 = x; 1622 for (nn=1; nn<p; nn++) { 1623 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1624 z0 = z1; 1625 z1 = z2; 1626 } 1627 Lpr = z2; 1628 A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(gllnodes[j]-gllnodes[r])*(gllnodes[j]-gllnodes[r])); 1629 } 1630 } 1631 } 1632 for (j=1; j<p+1; j++) { 1633 x = gllnodes[j]; 1634 z0 = 1.; 1635 z1 = x; 1636 for (nn=1; nn<p; nn++) { 1637 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1638 z0 = z1; 1639 z1 = z2; 1640 } 1641 Lpj = z2; 1642 A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+gllnodes[j])*(1.+gllnodes[j])); 1643 A[0][j] = A[j][0]; 1644 } 1645 for (j=0; j<p; j++) { 1646 x = gllnodes[j]; 1647 z0 = 1.; 1648 z1 = x; 1649 for (nn=1; nn<p; nn++) { 1650 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.)); 1651 z0 = z1; 1652 z1 = z2; 1653 } 1654 Lpj=z2; 1655 1656 A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-gllnodes[j])*(1.-gllnodes[j])); 1657 A[j][p] = A[p][j]; 1658 } 1659 A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.; 1660 A[p][p]=A[0][0]; 1661 *AA = A; 1662 PetscFunctionReturn(0); 1663 } 1664 1665 /*@C 1666 PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element 1667 1668 Not Collective 1669 1670 Input Parameter: 1671 + n - the number of GLL nodes 1672 . nodes - the GLL nodes 1673 . weights - the GLL weightss 1674 - A - the stiffness element 1675 1676 Level: beginner 1677 1678 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate() 1679 1680 @*/ 1681 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1682 { 1683 PetscErrorCode ierr; 1684 1685 PetscFunctionBegin; 1686 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 1687 ierr = PetscFree(*AA);CHKERRQ(ierr); 1688 *AA = NULL; 1689 PetscFunctionReturn(0); 1690 } 1691 1692 /*@C 1693 PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element 1694 1695 Not Collective 1696 1697 Input Parameter: 1698 + n - the number of GLL nodes 1699 . nodes - the GLL nodes 1700 . weights - the GLL weights 1701 1702 Output Parameter: 1703 . AA - the stiffness element 1704 - AAT - the transpose of AA (pass in NULL if you do not need this array) 1705 1706 Level: beginner 1707 1708 Notes: 1709 Destroy this with PetscGaussLobattoLegendreElementGradientDestroy() 1710 1711 You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented 1712 1713 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy() 1714 1715 @*/ 1716 PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT) 1717 { 1718 PetscReal **A, **AT = NULL; 1719 PetscErrorCode ierr; 1720 const PetscReal *gllnodes = nodes; 1721 const PetscInt p = n-1; 1722 PetscReal q,qp,Li, Lj,d0; 1723 PetscInt i,j; 1724 1725 PetscFunctionBegin; 1726 ierr = PetscMalloc1(n,&A);CHKERRQ(ierr); 1727 ierr = PetscMalloc1(n*n,&A[0]);CHKERRQ(ierr); 1728 for (i=1; i<n; i++) A[i] = A[i-1]+n; 1729 1730 if (AAT) { 1731 ierr = PetscMalloc1(n,&AT);CHKERRQ(ierr); 1732 ierr = PetscMalloc1(n*n,&AT[0]);CHKERRQ(ierr); 1733 for (i=1; i<n; i++) AT[i] = AT[i-1]+n; 1734 } 1735 1736 if (n==1) {A[0][0] = 0.;} 1737 d0 = (PetscReal)p*((PetscReal)p+1.)/4.; 1738 for (i=0; i<n; i++) { 1739 for (j=0; j<n; j++) { 1740 A[i][j] = 0.; 1741 qAndLEvaluation(p,gllnodes[i],&q,&qp,&Li); 1742 qAndLEvaluation(p,gllnodes[j],&q,&qp,&Lj); 1743 if (i!=j) A[i][j] = Li/(Lj*(gllnodes[i]-gllnodes[j])); 1744 if ((j==i) && (i==0)) A[i][j] = -d0; 1745 if (j==i && i==p) A[i][j] = d0; 1746 if (AT) AT[j][i] = A[i][j]; 1747 } 1748 } 1749 if (AAT) *AAT = AT; 1750 *AA = A; 1751 PetscFunctionReturn(0); 1752 } 1753 1754 /*@C 1755 PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGaussLobattoLegendreElementGradientCreate() 1756 1757 Not Collective 1758 1759 Input Parameter: 1760 + n - the number of GLL nodes 1761 . nodes - the GLL nodes 1762 . weights - the GLL weights 1763 . AA - the stiffness element 1764 - AAT - the transpose of the element 1765 1766 Level: beginner 1767 1768 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionCreate() 1769 1770 @*/ 1771 PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT) 1772 { 1773 PetscErrorCode ierr; 1774 1775 PetscFunctionBegin; 1776 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 1777 ierr = PetscFree(*AA);CHKERRQ(ierr); 1778 *AA = NULL; 1779 if (*AAT) { 1780 ierr = PetscFree((*AAT)[0]);CHKERRQ(ierr); 1781 ierr = PetscFree(*AAT);CHKERRQ(ierr); 1782 *AAT = NULL; 1783 } 1784 PetscFunctionReturn(0); 1785 } 1786 1787 /*@C 1788 PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element 1789 1790 Not Collective 1791 1792 Input Parameter: 1793 + n - the number of GLL nodes 1794 . nodes - the GLL nodes 1795 . weights - the GLL weightss 1796 1797 Output Parameter: 1798 . AA - the stiffness element 1799 1800 Level: beginner 1801 1802 Notes: 1803 Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy() 1804 1805 This is the same as the Gradient operator multiplied by the diagonal mass matrix 1806 1807 You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented 1808 1809 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionDestroy() 1810 1811 @*/ 1812 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1813 { 1814 PetscReal **D; 1815 PetscErrorCode ierr; 1816 const PetscReal *gllweights = weights; 1817 const PetscInt glln = n; 1818 PetscInt i,j; 1819 1820 PetscFunctionBegin; 1821 ierr = PetscGaussLobattoLegendreElementGradientCreate(n,nodes,weights,&D,NULL);CHKERRQ(ierr); 1822 for (i=0; i<glln; i++){ 1823 for (j=0; j<glln; j++) { 1824 D[i][j] = gllweights[i]*D[i][j]; 1825 } 1826 } 1827 *AA = D; 1828 PetscFunctionReturn(0); 1829 } 1830 1831 /*@C 1832 PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element 1833 1834 Not Collective 1835 1836 Input Parameter: 1837 + n - the number of GLL nodes 1838 . nodes - the GLL nodes 1839 . weights - the GLL weights 1840 - A - advection 1841 1842 Level: beginner 1843 1844 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementAdvectionCreate() 1845 1846 @*/ 1847 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1848 { 1849 PetscErrorCode ierr; 1850 1851 PetscFunctionBegin; 1852 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 1853 ierr = PetscFree(*AA);CHKERRQ(ierr); 1854 *AA = NULL; 1855 PetscFunctionReturn(0); 1856 } 1857 1858 PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1859 { 1860 PetscReal **A; 1861 PetscErrorCode ierr; 1862 const PetscReal *gllweights = weights; 1863 const PetscInt glln = n; 1864 PetscInt i,j; 1865 1866 PetscFunctionBegin; 1867 ierr = PetscMalloc1(glln,&A);CHKERRQ(ierr); 1868 ierr = PetscMalloc1(glln*glln,&A[0]);CHKERRQ(ierr); 1869 for (i=1; i<glln; i++) A[i] = A[i-1]+glln; 1870 if (glln==1) {A[0][0] = 0.;} 1871 for (i=0; i<glln; i++) { 1872 for (j=0; j<glln; j++) { 1873 A[i][j] = 0.; 1874 if (j==i) A[i][j] = gllweights[i]; 1875 } 1876 } 1877 *AA = A; 1878 PetscFunctionReturn(0); 1879 } 1880 1881 PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA) 1882 { 1883 PetscErrorCode ierr; 1884 1885 PetscFunctionBegin; 1886 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr); 1887 ierr = PetscFree(*AA);CHKERRQ(ierr); 1888 *AA = NULL; 1889 PetscFunctionReturn(0); 1890 } 1891