Actual source code: dt.c

petsc-3.13.3 2020-07-01
Report Typos and Errors
  1: /* Discretization tools */

  3:  #include <petscdt.h>
  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>

 11: #if defined(PETSC_HAVE_MPFR)
 12: #include <mpfr.h>
 13: #endif

 15: const char *const PetscDTNodeTypes[] = {"gaussjacobi", "equispaced", "tanhsinh", "PETSCDTNODES_", 0};

 17: static PetscBool GolubWelschCite       = PETSC_FALSE;
 18: const char       GolubWelschCitation[] = "@article{GolubWelsch1969,\n"
 19:                                          "  author  = {Golub and Welsch},\n"
 20:                                          "  title   = {Calculation of Quadrature Rules},\n"
 21:                                          "  journal = {Math. Comp.},\n"
 22:                                          "  volume  = {23},\n"
 23:                                          "  number  = {106},\n"
 24:                                          "  pages   = {221--230},\n"
 25:                                          "  year    = {1969}\n}\n";

 27: /* Numerical tests in src/dm/dt/tests/ex1.c show that when computing the nodes and weights of Gauss-Jacobi
 28:    quadrature rules:

 30:    - in double precision, Newton's method and Golub & Welsch both work for moderate degrees (< 100),
 31:    - in single precision, Newton's method starts producing incorrect roots around n = 15, but
 32:      the weights from Golub & Welsch become a problem before then: they produces errors
 33:      in computing the Jacobi-polynomial Gram matrix around n = 6.

 35:    So we default to Newton's method (required fewer dependencies) */
 36: PetscBool PetscDTGaussQuadratureNewton_Internal = PETSC_TRUE;

 38: PetscClassId PETSCQUADRATURE_CLASSID = 0;

 40: /*@
 41:   PetscQuadratureCreate - Create a PetscQuadrature object

 43:   Collective

 45:   Input Parameter:
 46: . comm - The communicator for the PetscQuadrature object

 48:   Output Parameter:
 49: . q  - The PetscQuadrature object

 51:   Level: beginner

 53: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
 54: @*/
 55: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
 56: {

 61:   DMInitializePackage();
 62:   PetscHeaderCreate(*q,PETSCQUADRATURE_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
 63:   (*q)->dim       = -1;
 64:   (*q)->Nc        =  1;
 65:   (*q)->order     = -1;
 66:   (*q)->numPoints = 0;
 67:   (*q)->points    = NULL;
 68:   (*q)->weights   = NULL;
 69:   return(0);
 70: }

 72: /*@
 73:   PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object

 75:   Collective on q

 77:   Input Parameter:
 78: . q  - The PetscQuadrature object

 80:   Output Parameter:
 81: . r  - The new PetscQuadrature object

 83:   Level: beginner

 85: .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData()
 86: @*/
 87: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
 88: {
 89:   PetscInt         order, dim, Nc, Nq;
 90:   const PetscReal *points, *weights;
 91:   PetscReal       *p, *w;
 92:   PetscErrorCode   ierr;

 96:   PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);
 97:   PetscQuadratureGetOrder(q, &order);
 98:   PetscQuadratureSetOrder(*r, order);
 99:   PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);
100:   PetscMalloc1(Nq*dim, &p);
101:   PetscMalloc1(Nq*Nc, &w);
102:   PetscArraycpy(p, points, Nq*dim);
103:   PetscArraycpy(w, weights, Nc * Nq);
104:   PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);
105:   return(0);
106: }

108: /*@
109:   PetscQuadratureDestroy - Destroys a PetscQuadrature object

111:   Collective on q

113:   Input Parameter:
114: . q  - The PetscQuadrature object

116:   Level: beginner

118: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
119: @*/
120: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
121: {

125:   if (!*q) return(0);
127:   if (--((PetscObject)(*q))->refct > 0) {
128:     *q = NULL;
129:     return(0);
130:   }
131:   PetscFree((*q)->points);
132:   PetscFree((*q)->weights);
133:   PetscHeaderDestroy(q);
134:   return(0);
135: }

137: /*@
138:   PetscQuadratureGetOrder - Return the order of the method

140:   Not collective

142:   Input Parameter:
143: . q - The PetscQuadrature object

145:   Output Parameter:
146: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated

148:   Level: intermediate

150: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
151: @*/
152: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
153: {
157:   *order = q->order;
158:   return(0);
159: }

161: /*@
162:   PetscQuadratureSetOrder - Return the order of the method

164:   Not collective

166:   Input Parameters:
167: + q - The PetscQuadrature object
168: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated

170:   Level: intermediate

172: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
173: @*/
174: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
175: {
178:   q->order = order;
179:   return(0);
180: }

182: /*@
183:   PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated

185:   Not collective

187:   Input Parameter:
188: . q - The PetscQuadrature object

190:   Output Parameter:
191: . Nc - The number of components

193:   Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.

195:   Level: intermediate

197: .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
198: @*/
199: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
200: {
204:   *Nc = q->Nc;
205:   return(0);
206: }

208: /*@
209:   PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated

211:   Not collective

213:   Input Parameters:
214: + q  - The PetscQuadrature object
215: - Nc - The number of components

217:   Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.

219:   Level: intermediate

221: .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
222: @*/
223: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
224: {
227:   q->Nc = Nc;
228:   return(0);
229: }

231: /*@C
232:   PetscQuadratureGetData - Returns the data defining the quadrature

234:   Not collective

236:   Input Parameter:
237: . q  - The PetscQuadrature object

239:   Output Parameters:
240: + dim - The spatial dimension
241: . Nc - The number of components
242: . npoints - The number of quadrature points
243: . points - The coordinates of each quadrature point
244: - weights - The weight of each quadrature point

246:   Level: intermediate

248:   Fortran Notes:
249:     From Fortran you must call PetscQuadratureRestoreData() when you are done with the data

251: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
252: @*/
253: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
254: {
257:   if (dim) {
259:     *dim = q->dim;
260:   }
261:   if (Nc) {
263:     *Nc = q->Nc;
264:   }
265:   if (npoints) {
267:     *npoints = q->numPoints;
268:   }
269:   if (points) {
271:     *points = q->points;
272:   }
273:   if (weights) {
275:     *weights = q->weights;
276:   }
277:   return(0);
278: }

280: static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[])
281: {
282:   PetscScalar    *Js, *Jinvs;
283:   PetscInt       i, j, k;
284:   PetscBLASInt   bm, bn, info;

288:   if (!m || !n) return(0);
289:   PetscBLASIntCast(m, &bm);
290:   PetscBLASIntCast(n, &bn);
291: #if defined(PETSC_USE_COMPLEX)
292:   PetscMalloc2(m*n, &Js, m*n, &Jinvs);
293:   for (i = 0; i < m*n; i++) Js[i] = J[i];
294: #else
295:   Js = (PetscReal *) J;
296:   Jinvs = Jinv;
297: #endif
298:   if (m == n) {
299:     PetscBLASInt *pivots;
300:     PetscScalar *W;

302:     PetscMalloc2(m, &pivots, m, &W);

304:     PetscArraycpy(Jinvs, Js, m * m);
305:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info));
306:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
307:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &info));
308:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
309:     PetscFree2(pivots, W);
310:   } else if (m < n) {
311:     PetscScalar *JJT;
312:     PetscBLASInt *pivots;
313:     PetscScalar *W;

315:     PetscMalloc1(m*m, &JJT);
316:     PetscMalloc2(m, &pivots, m, &W);
317:     for (i = 0; i < m; i++) {
318:       for (j = 0; j < m; j++) {
319:         PetscScalar val = 0.;

321:         for (k = 0; k < n; k++) val += Js[i * n + k] * Js[j * n + k];
322:         JJT[i * m + j] = val;
323:       }
324:     }

326:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info));
327:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
328:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info));
329:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
330:     for (i = 0; i < n; i++) {
331:       for (j = 0; j < m; j++) {
332:         PetscScalar val = 0.;

334:         for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j];
335:         Jinvs[i * m + j] = val;
336:       }
337:     }
338:     PetscFree2(pivots, W);
339:     PetscFree(JJT);
340:   } else {
341:     PetscScalar *JTJ;
342:     PetscBLASInt *pivots;
343:     PetscScalar *W;

345:     PetscMalloc1(n*n, &JTJ);
346:     PetscMalloc2(n, &pivots, n, &W);
347:     for (i = 0; i < n; i++) {
348:       for (j = 0; j < n; j++) {
349:         PetscScalar val = 0.;

351:         for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j];
352:         JTJ[i * n + j] = val;
353:       }
354:     }

356:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bn, pivots, &info));
357:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
358:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info));
359:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
360:     for (i = 0; i < n; i++) {
361:       for (j = 0; j < m; j++) {
362:         PetscScalar val = 0.;

364:         for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k];
365:         Jinvs[i * m + j] = val;
366:       }
367:     }
368:     PetscFree2(pivots, W);
369:     PetscFree(JTJ);
370:   }
371: #if defined(PETSC_USE_COMPLEX)
372:   for (i = 0; i < m*n; i++) Jinv[i] = PetscRealPart(Jinvs[i]);
373:   PetscFree2(Js, Jinvs);
374: #endif
375:   return(0);
376: }

378: /*@
379:    PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation.

381:    Collecive on PetscQuadrature

383:    Input Arguments:
384: +  q - the quadrature functional
385: .  imageDim - the dimension of the image of the transformation
386: .  origin - a point in the original space
387: .  originImage - the image of the origin under the transformation
388: .  J - the Jacobian of the image: an [imageDim x dim] matrix in row major order
389: -  formDegree - transform the quadrature weights as k-forms of this form degree (if the number of components is a multiple of (dim choose formDegree), it is assumed that they represent multiple k-forms) [see PetscDTAltVPullback() for interpretation of formDegree]

391:    Output Arguments:
392: .  Jinvstarq - a quadrature rule where each point is the image of a point in the original quadrature rule, and where the k-form weights have been pulled-back by the pseudoinverse of J to the k-form weights in the image space.

394:    Note: the new quadrature rule will have a different number of components if spaces have different dimensions.  For example, pushing a 2-form forward from a two dimensional space to a three dimensional space changes the number of components from 1 to 3.

396:    Level: intermediate

398: .seealso: PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
399: @*/
400: PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq)
401: {
402:   PetscInt         dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c;
403:   const PetscReal *points;
404:   const PetscReal *weights;
405:   PetscReal       *imagePoints, *imageWeights;
406:   PetscReal       *Jinv;
407:   PetscReal       *Jinvstar;
408:   PetscErrorCode   ierr;

412:   if (imageDim < PetscAbsInt(formDegree)) SETERRQ2(PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Cannot represent a %D-form in %D dimensions", PetscAbsInt(formDegree), imageDim);
413:   PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights);
414:   PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize);
415:   if (Nc % formSize) SETERRQ2(PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Number of components %D is not a multiple of formSize %D\n", Nc, formSize);
416:   Ncopies = Nc / formSize;
417:   PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize);
418:   imageNc = Ncopies * imageFormSize;
419:   PetscMalloc1(Npoints * imageDim, &imagePoints);
420:   PetscMalloc1(Npoints * imageNc, &imageWeights);
421:   PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar);
422:   PetscDTJacobianInverse_Internal(imageDim, dim, J, Jinv);
423:   PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar);
424:   for (pt = 0; pt < Npoints; pt++) {
425:     const PetscReal *point = &points[pt * dim];
426:     PetscReal       *imagePoint = &imagePoints[pt * imageDim];

428:     for (i = 0; i < imageDim; i++) {
429:       PetscReal val = originImage[i];

431:       for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]);
432:       imagePoint[i] = val;
433:     }
434:     for (c = 0; c < Ncopies; c++) {
435:       const PetscReal *form = &weights[pt * Nc + c * formSize];
436:       PetscReal       *imageForm = &imageWeights[pt * imageNc + c * imageFormSize];

438:       for (i = 0; i < imageFormSize; i++) {
439:         PetscReal val = 0.;

441:         for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j];
442:         imageForm[i] = val;
443:       }
444:     }
445:   }
446:   PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq);
447:   PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights);
448:   PetscFree2(Jinv, Jinvstar);
449:   return(0);
450: }

452: /*@C
453:   PetscQuadratureSetData - Sets the data defining the quadrature

455:   Not collective

457:   Input Parameters:
458: + q  - The PetscQuadrature object
459: . dim - The spatial dimension
460: . Nc - The number of components
461: . npoints - The number of quadrature points
462: . points - The coordinates of each quadrature point
463: - weights - The weight of each quadrature point

465:   Note: This routine owns the references to points and weights, so they must be allocated using PetscMalloc() and the user should not free them.

467:   Level: intermediate

469: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
470: @*/
471: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
472: {
475:   if (dim >= 0)     q->dim       = dim;
476:   if (Nc >= 0)      q->Nc        = Nc;
477:   if (npoints >= 0) q->numPoints = npoints;
478:   if (points) {
480:     q->points = points;
481:   }
482:   if (weights) {
484:     q->weights = weights;
485:   }
486:   return(0);
487: }

489: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
490: {
491:   PetscInt          q, d, c;
492:   PetscViewerFormat format;
493:   PetscErrorCode    ierr;

496:   if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D) with %D components\n", quad->order, quad->numPoints, quad->dim, quad->Nc);}
497:   else              {PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D)\n", quad->order, quad->numPoints, quad->dim);}
498:   PetscViewerGetFormat(v, &format);
499:   if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
500:   for (q = 0; q < quad->numPoints; ++q) {
501:     PetscViewerASCIIPrintf(v, "p%D (", q);
502:     PetscViewerASCIIUseTabs(v, PETSC_FALSE);
503:     for (d = 0; d < quad->dim; ++d) {
504:       if (d) {PetscViewerASCIIPrintf(v, ", ");}
505:       PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q*quad->dim+d]);
506:     }
507:     PetscViewerASCIIPrintf(v, ") ");
508:     if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "w%D (", q);}
509:     for (c = 0; c < quad->Nc; ++c) {
510:       if (c) {PetscViewerASCIIPrintf(v, ", ");}
511:       PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q*quad->Nc+c]);
512:     }
513:     if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, ")");}
514:     PetscViewerASCIIPrintf(v, "\n");
515:     PetscViewerASCIIUseTabs(v, PETSC_TRUE);
516:   }
517:   return(0);
518: }

520: /*@C
521:   PetscQuadratureView - Views a PetscQuadrature object

523:   Collective on quad

525:   Input Parameters:
526: + quad  - The PetscQuadrature object
527: - viewer - The PetscViewer object

529:   Level: beginner

531: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
532: @*/
533: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
534: {
535:   PetscBool      iascii;

541:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) quad), &viewer);}
542:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
543:   PetscViewerASCIIPushTab(viewer);
544:   if (iascii) {PetscQuadratureView_Ascii(quad, viewer);}
545:   PetscViewerASCIIPopTab(viewer);
546:   return(0);
547: }

549: /*@C
550:   PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement

552:   Not collective

554:   Input Parameter:
555: + q - The original PetscQuadrature
556: . numSubelements - The number of subelements the original element is divided into
557: . v0 - An array of the initial points for each subelement
558: - jac - An array of the Jacobian mappings from the reference to each subelement

560:   Output Parameters:
561: . dim - The dimension

563:   Note: Together v0 and jac define an affine mapping from the original reference element to each subelement

565:  Not available from Fortran

567:   Level: intermediate

569: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
570: @*/
571: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
572: {
573:   const PetscReal *points,    *weights;
574:   PetscReal       *pointsRef, *weightsRef;
575:   PetscInt         dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
576:   PetscErrorCode   ierr;

583:   PetscQuadratureCreate(PETSC_COMM_SELF, qref);
584:   PetscQuadratureGetOrder(q, &order);
585:   PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
586:   npointsRef = npoints*numSubelements;
587:   PetscMalloc1(npointsRef*dim,&pointsRef);
588:   PetscMalloc1(npointsRef*Nc, &weightsRef);
589:   for (c = 0; c < numSubelements; ++c) {
590:     for (p = 0; p < npoints; ++p) {
591:       for (d = 0; d < dim; ++d) {
592:         pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
593:         for (e = 0; e < dim; ++e) {
594:           pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
595:         }
596:       }
597:       /* Could also use detJ here */
598:       for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
599:     }
600:   }
601:   PetscQuadratureSetOrder(*qref, order);
602:   PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
603:   return(0);
604: }

606: /* Compute the coefficients for the Jacobi polynomial recurrence,
607:  *
608:  * J^{a,b}_n(x) = (cnm1 + cnm1x * x) * J^{a,b}_{n-1}(x) - cnm2 * J^{a,b}_{n-2}(x).
609:  */
610: #define PetscDTJacobiRecurrence_Internal(n,a,b,cnm1,cnm1x,cnm2) \
611: do {                                                            \
612:   PetscReal _a = (a);                                           \
613:   PetscReal _b = (b);                                           \
614:   PetscReal _n = (n);                                           \
615:   if (n == 1) {                                                 \
616:     (cnm1) = (_a-_b) * 0.5;                                     \
617:     (cnm1x) = (_a+_b+2.)*0.5;                                   \
618:     (cnm2) = 0.;                                                \
619:   } else {                                                      \
620:     PetscReal _2n = _n+_n;                                      \
621:     PetscReal _d = (_2n*(_n+_a+_b)*(_2n+_a+_b-2));              \
622:     PetscReal _n1 = (_2n+_a+_b-1.)*(_a*_a-_b*_b);               \
623:     PetscReal _n1x = (_2n+_a+_b-1.)*(_2n+_a+_b)*(_2n+_a+_b-2);  \
624:     PetscReal _n2 = 2.*((_n+_a-1.)*(_n+_b-1.)*(_2n+_a+_b));     \
625:     (cnm1) = _n1 / _d;                                          \
626:     (cnm1x) = _n1x / _d;                                        \
627:     (cnm2) = _n2 / _d;                                          \
628:   }                                                             \
629: } while (0)

631: static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p)
632: {
633:   PetscReal ak, bk;
634:   PetscReal abk1;
635:   PetscInt i,l,maxdegree;

638:   maxdegree = degrees[ndegree-1] - k;
639:   ak = a + k;
640:   bk = b + k;
641:   abk1 = a + b + k + 1.;
642:   if (maxdegree < 0) {
643:     for (i = 0; i < npoints; i++) for (l = 0; l < ndegree; l++) p[i*ndegree+l] = 0.;
644:     return(0);
645:   }
646:   for (i=0; i<npoints; i++) {
647:     PetscReal pm1,pm2,x;
648:     PetscReal cnm1, cnm1x, cnm2;
649:     PetscInt  j,m;

651:     x    = points[i];
652:     pm2  = 1.;
653:     PetscDTJacobiRecurrence_Internal(1,ak,bk,cnm1,cnm1x,cnm2);
654:     pm1 = (cnm1 + cnm1x*x);
655:     l    = 0;
656:     while (l < ndegree && degrees[l] - k < 0) {
657:       p[l++] = 0.;
658:     }
659:     while (l < ndegree && degrees[l] - k == 0) {
660:       p[l] = pm2;
661:       for (m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5;
662:       l++;
663:     }
664:     while (l < ndegree && degrees[l] - k == 1) {
665:       p[l] = pm1;
666:       for (m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5;
667:       l++;
668:     }
669:     for (j=2; j<=maxdegree; j++) {
670:       PetscReal pp;

672:       PetscDTJacobiRecurrence_Internal(j,ak,bk,cnm1,cnm1x,cnm2);
673:       pp   = (cnm1 + cnm1x*x)*pm1 - cnm2*pm2;
674:       pm2  = pm1;
675:       pm1  = pp;
676:       while (l < ndegree && degrees[l] - k == j) {
677:         p[l] = pp;
678:         for (m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5;
679:         l++;
680:       }
681:     }
682:     p += ndegree;
683:   }
684:   return(0);
685: }

687: /*@
688:    PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$
689:                        at points

691:    Not Collective

693:    Input Arguments:
694: +  npoints - number of spatial points to evaluate at
695: .  alpha - the left exponent > -1
696: .  beta - the right exponent > -1
697: .  points - array of locations to evaluate at
698: .  ndegree - number of basis degrees to evaluate
699: -  degrees - sorted array of degrees to evaluate

701:    Output Arguments:
702: +  B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
703: .  D - row-oriented derivative evaluation matrix (or NULL)
704: -  D2 - row-oriented second derivative evaluation matrix (or NULL)

706:    Level: intermediate

708: .seealso: PetscDTGaussQuadrature()
709: @*/
710: PetscErrorCode PetscDTJacobiEval(PetscInt npoints,PetscReal alpha, PetscReal beta, const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
711: {

715:   if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1.");
716:   if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1.");
717:   if (!npoints || !ndegree) return(0);
718:   if (B)  {PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B);}
719:   if (D)  {PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D);}
720:   if (D2) {PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2);}
721:   return(0);
722: }

724: /*@
725:    PetscDTLegendreEval - evaluate Legendre polynomials at points

727:    Not Collective

729:    Input Arguments:
730: +  npoints - number of spatial points to evaluate at
731: .  points - array of locations to evaluate at
732: .  ndegree - number of basis degrees to evaluate
733: -  degrees - sorted array of degrees to evaluate

735:    Output Arguments:
736: +  B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
737: .  D - row-oriented derivative evaluation matrix (or NULL)
738: -  D2 - row-oriented second derivative evaluation matrix (or NULL)

740:    Level: intermediate

742: .seealso: PetscDTGaussQuadrature()
743: @*/
744: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
745: {

749:   PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2);
750:   return(0);
751: }

753: /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V
754:  * with lds n; diag and subdiag are overwritten */
755: static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[],
756:                                                             PetscReal eigs[], PetscScalar V[])
757: {
758:   char jobz = 'V'; /* eigenvalues and eigenvectors */
759:   char range = 'A'; /* all eigenvalues will be found */
760:   PetscReal VL = 0.; /* ignored because range is 'A' */
761:   PetscReal VU = 0.; /* ignored because range is 'A' */
762:   PetscBLASInt IL = 0; /* ignored because range is 'A' */
763:   PetscBLASInt IU = 0; /* ignored because range is 'A' */
764:   PetscReal abstol = 0.; /* unused */
765:   PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */
766:   PetscBLASInt *isuppz;
767:   PetscBLASInt lwork, liwork;
768:   PetscReal workquery;
769:   PetscBLASInt  iworkquery;
770:   PetscBLASInt *iwork;
771:   PetscBLASInt info;
772:   PetscReal *work = NULL;

776: #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG)
777:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
778: #endif
779:   PetscBLASIntCast(n, &bn);
780:   PetscBLASIntCast(n, &ldz);
781: #if !defined(PETSC_MISSING_LAPACK_STEGR)
782:   PetscMalloc1(2 * n, &isuppz);
783:   lwork = -1;
784:   liwork = -1;
785:   PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,&workquery,&lwork,&iworkquery,&liwork,&info));
786:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEGR error");
787:   lwork = (PetscBLASInt) workquery;
788:   liwork = (PetscBLASInt) iworkquery;
789:   PetscMalloc2(lwork, &work, liwork, &iwork);
790:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
791:   PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,work,&lwork,iwork,&liwork,&info));
792:   PetscFPTrapPop();
793:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEGR error");
794:   PetscFree2(work, iwork);
795:   PetscFree(isuppz);
796: #elif !defined(PETSC_MISSING_LAPACK_STEQR)
797:   jobz = 'I'; /* Compute eigenvalues and eigenvectors of the
798:                  tridiagonal matrix.  Z is initialized to the identity
799:                  matrix. */
800:   PetscMalloc1(PetscMax(1,2*n-2),&work);
801:   PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&bn,diag,subdiag,V,&ldz,work,&info));
802:   PetscFPTrapPop();
803:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");
804:   PetscFree(work);
805:   PetscArraycpy(eigs,diag,n);
806: #endif
807:   return(0);
808: }

810: /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi
811:  * quadrature rules on the interval [-1, 1] */
812: static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw)
813: {
814:   PetscReal twoab1;
815:   PetscInt  m = n - 2;
816:   PetscReal a = alpha + 1.;
817:   PetscReal b = beta + 1.;
818:   PetscReal gra, grb;

821:   twoab1 = PetscPowReal(2., a + b - 1.);
822: #if defined(PETSC_HAVE_LGAMMA)
823:   grb = PetscExpReal(2. * PetscLGamma(b+1.) + PetscLGamma(m+1.) + PetscLGamma(m+a+1.) -
824:                      (PetscLGamma(m+b+1) + PetscLGamma(m+a+b+1.)));
825:   gra = PetscExpReal(2. * PetscLGamma(a+1.) + PetscLGamma(m+1.) + PetscLGamma(m+b+1.) -
826:                      (PetscLGamma(m+a+1) + PetscLGamma(m+a+b+1.)));
827: #else
828:   {
829:     PetscInt alphai = (PetscInt) alpha;
830:     PetscInt betai = (PetscInt) beta;

833:     if ((PetscReal) alphai == alpha && (PetscReal) betai == beta) {
834:       PetscReal binom1, binom2;

836:       PetscDTBinomial(m+b, b, &binom1);
837:       PetscDTBinomial(m+a+b, b, &binom2);
838:       grb = 1./ (binom1 * binom2);
839:       PetscDTBinomial(m+a, a, &binom1);
840:       PetscDTBinomial(m+a+b, a, &binom2);
841:       gra = 1./ (binom1 * binom2);
842:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
843:   }
844: #endif
845:   *leftw = twoab1 * grb / b;
846:   *rightw = twoab1 * gra / a;
847:   return(0);
848: }

850: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
851:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
852: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
853: {
854:   PetscReal pn1, pn2;
855:   PetscReal cnm1, cnm1x, cnm2;
856:   PetscInt  k;

859:   if (!n) {*P = 1.0; return(0);}
860:   PetscDTJacobiRecurrence_Internal(1,a,b,cnm1,cnm1x,cnm2);
861:   pn2 = 1.;
862:   pn1 = cnm1 + cnm1x*x;
863:   if (n == 1) {*P = pn1; return(0);}
864:   *P  = 0.0;
865:   for (k = 2; k < n+1; ++k) {
866:     PetscDTJacobiRecurrence_Internal(k,a,b,cnm1,cnm1x,cnm2);

868:     *P  = (cnm1 + cnm1x*x)*pn1 - cnm2*pn2;
869:     pn2 = pn1;
870:     pn1 = *P;
871:   }
872:   return(0);
873: }

875: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
876: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P)
877: {
878:   PetscReal      nP;
879:   PetscInt       i;

883:   *P = 0.0;
884:   if (k > n) return(0);
885:   PetscDTComputeJacobi(a+k, b+k, n-k, x, &nP);
886:   for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5;
887:   *P = nP;
888:   return(0);
889: }

891: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
892: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
893: {
895:   *xi  = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
896:   *eta = y;
897:   return(0);
898: }

900: static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[])
901: {
902:   PetscInt       maxIter = 100;
903:   PetscReal      eps     = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON));
904:   PetscReal      a1, a6, gf;
905:   PetscInt       k;


910:   a1 = PetscPowReal(2.0, a+b+1);
911: #if defined(PETSC_HAVE_LGAMMA)
912:   {
913:     PetscReal a2, a3, a4, a5;
914:     a2 = PetscLGamma(a + npoints + 1);
915:     a3 = PetscLGamma(b + npoints + 1);
916:     a4 = PetscLGamma(a + b + npoints + 1);
917:     a5 = PetscLGamma(npoints + 1);
918:     gf = PetscExpReal(a2 + a3 - (a4 + a5));
919:   }
920: #else
921:   {
922:     PetscInt ia, ib;

924:     ia = (PetscInt) a;
925:     ib = (PetscInt) b;
926:     gf = 1.;
927:     if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */
928:       for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k);
929:     } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */
930:       for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k);
931:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
932:   }
933: #endif

935:   a6   = a1 * gf;
936:   /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
937:    Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
938:   for (k = 0; k < npoints; ++k) {
939:     PetscReal r = PetscCosReal(PETSC_PI * (1. - (4.*k + 3. + 2.*b) / (4.*npoints + 2.*(a + b + 1.)))), dP;
940:     PetscInt  j;

942:     if (k > 0) r = 0.5 * (r + x[k-1]);
943:     for (j = 0; j < maxIter; ++j) {
944:       PetscReal s = 0.0, delta, f, fp;
945:       PetscInt  i;

947:       for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
948:       PetscDTComputeJacobi(a, b, npoints, r, &f);
949:       PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp);
950:       delta = f / (fp - f * s);
951:       r     = r - delta;
952:       if (PetscAbsReal(delta) < eps) break;
953:     }
954:     x[k] = r;
955:     PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP);
956:     w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
957:   }
958:   return(0);
959: }

961: /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi
962:  * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */
963: static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s)
964: {
965:   PetscInt       i;

968:   for (i = 0; i < nPoints; i++) {
969:     PetscReal A, B, C;

971:     PetscDTJacobiRecurrence_Internal(i+1,a,b,A,B,C);
972:     d[i] = -A / B;
973:     if (i) s[i-1] *= C / B;
974:     if (i < nPoints - 1) s[i] = 1. / B;
975:   }
976:   return(0);
977: }

979: static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
980: {
981:   PetscReal mu0;
982:   PetscReal ga, gb, gab;
983:   PetscInt i;

987:   PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite);

989: #if defined(PETSC_HAVE_TGAMMA)
990:   ga  = PetscTGamma(a + 1);
991:   gb  = PetscTGamma(b + 1);
992:   gab = PetscTGamma(a + b + 2);
993: #else
994:   {
995:     PetscInt ia, ib;

997:     ia = (PetscInt) a;
998:     ib = (PetscInt) b;
999:     if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */
1000:       PetscDTFactorial(ia, &ga);
1001:       PetscDTFactorial(ib, &gb);
1002:       PetscDTFactorial(ia + ib + 1, &gb);
1003:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
1004:   }
1005: #endif
1006:   mu0 = PetscPowReal(2.,a + b + 1.) * ga * gb / gab;

1008: #if defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1009:   {
1010:     PetscReal *diag, *subdiag;
1011:     PetscScalar *V;

1013:     PetscMalloc2(npoints, &diag, npoints, &subdiag);
1014:     PetscMalloc1(npoints*npoints, &V);
1015:     PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag);
1016:     for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]);
1017:     PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V);
1018:     for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0;
1019:     PetscFree(V);
1020:     PetscFree2(diag, subdiag);
1021:   }
1022: #else
1023:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1024: #endif
1025:   { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the
1026:        eigenvalues are not guaranteed to be in ascending order.  So we heave a passive aggressive sigh and check that
1027:        the eigenvalues are sorted */
1028:     PetscBool sorted;

1030:     PetscSortedReal(npoints, x, &sorted);
1031:     if (!sorted) {
1032:       PetscInt *order, i;
1033:       PetscReal *tmp;

1035:       PetscMalloc2(npoints, &order, npoints, &tmp);
1036:       for (i = 0; i < npoints; i++) order[i] = i;
1037:       PetscSortRealWithPermutation(npoints, x, order);
1038:       PetscArraycpy(tmp, x, npoints);
1039:       for (i = 0; i < npoints; i++) x[i] = tmp[order[i]];
1040:       PetscArraycpy(tmp, w, npoints);
1041:       for (i = 0; i < npoints; i++) w[i] = tmp[order[i]];
1042:       PetscFree2(order, tmp);
1043:     }
1044:   }
1045:   return(0);
1046: }

1048: static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1049: {

1053:   if (npoints < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be positive");
1054:   /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1055:   if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1.");
1056:   if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1.");

1058:   if (newton) {
1059:     PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w);
1060:   } else {
1061:     PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w);
1062:   }
1063:   if (alpha == beta) { /* symmetrize */
1064:     PetscInt i;
1065:     for (i = 0; i < (npoints + 1) / 2; i++) {
1066:       PetscInt  j  = npoints - 1 - i;
1067:       PetscReal xi = x[i];
1068:       PetscReal xj = x[j];
1069:       PetscReal wi = w[i];
1070:       PetscReal wj = w[j];

1072:       x[i] = (xi - xj) / 2.;
1073:       x[j] = (xj - xi) / 2.;
1074:       w[i] = w[j] = (wi + wj) / 2.;
1075:     }
1076:   }
1077:   return(0);
1078: }

1080: /*@
1081:   PetscDTGaussJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1082:   $(x-a)^\alpha (x-b)^\beta$.

1084:   Not collective

1086:   Input Parameters:
1087: + npoints - the number of points in the quadrature rule
1088: . a - the left endpoint of the interval
1089: . b - the right endpoint of the interval
1090: . alpha - the left exponent
1091: - beta - the right exponent

1093:   Output Parameters:
1094: + x - array of length npoints, the locations of the quadrature points
1095: - w - array of length npoints, the weights of the quadrature points

1097:   Level: intermediate

1099:   Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 1.
1100: @*/
1101: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1102: {
1103:   PetscInt       i;

1107:   PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);
1108:   if (a != -1. || b != 1.) { /* shift */
1109:     for (i = 0; i < npoints; i++) {
1110:       x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1111:       w[i] *= (b - a) / 2.;
1112:     }
1113:   }
1114:   return(0);
1115: }

1117: static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1118: {
1119:   PetscInt       i;

1123:   if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be positive");
1124:   /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1125:   if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1.");
1126:   if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1.");

1128:   x[0] = -1.;
1129:   x[npoints-1] = 1.;
1130:   if (npoints > 2) {
1131:     PetscDTGaussJacobiQuadrature_Internal(npoints-2, alpha+1., beta+1., &x[1], &w[1], newton);
1132:   }
1133:   for (i = 1; i < npoints - 1; i++) {
1134:     w[i] /= (1. - x[i]*x[i]);
1135:   }
1136:   PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints-1]);
1137:   return(0);
1138: }

1140: /*@
1141:   PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1142:   $(x-a)^\alpha (x-b)^\beta$, with endpoints a and b included as quadrature points.

1144:   Not collective

1146:   Input Parameters:
1147: + npoints - the number of points in the quadrature rule
1148: . a - the left endpoint of the interval
1149: . b - the right endpoint of the interval
1150: . alpha - the left exponent
1151: - beta - the right exponent

1153:   Output Parameters:
1154: + x - array of length npoints, the locations of the quadrature points
1155: - w - array of length npoints, the weights of the quadrature points

1157:   Level: intermediate

1159:   Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 3.
1160: @*/
1161: PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1162: {
1163:   PetscInt       i;

1167:   PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);
1168:   if (a != -1. || b != 1.) { /* shift */
1169:     for (i = 0; i < npoints; i++) {
1170:       x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1171:       w[i] *= (b - a) / 2.;
1172:     }
1173:   }
1174:   return(0);
1175: }

1177: /*@
1178:    PetscDTGaussQuadrature - create Gauss-Legendre quadrature

1180:    Not Collective

1182:    Input Arguments:
1183: +  npoints - number of points
1184: .  a - left end of interval (often-1)
1185: -  b - right end of interval (often +1)

1187:    Output Arguments:
1188: +  x - quadrature points
1189: -  w - quadrature weights

1191:    Level: intermediate

1193:    References:
1194: .   1. - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.

1196: .seealso: PetscDTLegendreEval()
1197: @*/
1198: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
1199: {
1200:   PetscInt       i;

1204:   PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal);
1205:   if (a != -1. || b != 1.) { /* shift */
1206:     for (i = 0; i < npoints; i++) {
1207:       x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1208:       w[i] *= (b - a) / 2.;
1209:     }
1210:   }
1211:   return(0);
1212: }

1214: /*@C
1215:    PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre
1216:                       nodes of a given size on the domain [-1,1]

1218:    Not Collective

1220:    Input Parameter:
1221: +  n - number of grid nodes
1222: -  type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON

1224:    Output Arguments:
1225: +  x - quadrature points
1226: -  w - quadrature weights

1228:    Notes:
1229:     For n > 30  the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
1230:           close enough to the desired solution

1232:    These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes

1234:    See  https://epubs.siam.org/doi/abs/10.1137/110855442  https://epubs.siam.org/doi/abs/10.1137/120889873 for better ways to compute GLL nodes

1236:    Level: intermediate

1238: .seealso: PetscDTGaussQuadrature()

1240: @*/
1241: PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal *x,PetscReal *w)
1242: {
1243:   PetscBool      newton;

1247:   if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide at least 2 grid points per element");
1248:   newton = (PetscBool) (type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON);
1249:   PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton);
1250:   return(0);
1251: }

1253: /*@
1254:   PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature

1256:   Not Collective

1258:   Input Arguments:
1259: + dim     - The spatial dimension
1260: . Nc      - The number of components
1261: . npoints - number of points in one dimension
1262: . a       - left end of interval (often-1)
1263: - b       - right end of interval (often +1)

1265:   Output Argument:
1266: . q - A PetscQuadrature object

1268:   Level: intermediate

1270: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
1271: @*/
1272: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1273: {
1274:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
1275:   PetscReal     *x, *w, *xw, *ww;

1279:   PetscMalloc1(totpoints*dim,&x);
1280:   PetscMalloc1(totpoints*Nc,&w);
1281:   /* Set up the Golub-Welsch system */
1282:   switch (dim) {
1283:   case 0:
1284:     PetscFree(x);
1285:     PetscFree(w);
1286:     PetscMalloc1(1, &x);
1287:     PetscMalloc1(Nc, &w);
1288:     x[0] = 0.0;
1289:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
1290:     break;
1291:   case 1:
1292:     PetscMalloc1(npoints,&ww);
1293:     PetscDTGaussQuadrature(npoints, a, b, x, ww);
1294:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
1295:     PetscFree(ww);
1296:     break;
1297:   case 2:
1298:     PetscMalloc2(npoints,&xw,npoints,&ww);
1299:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
1300:     for (i = 0; i < npoints; ++i) {
1301:       for (j = 0; j < npoints; ++j) {
1302:         x[(i*npoints+j)*dim+0] = xw[i];
1303:         x[(i*npoints+j)*dim+1] = xw[j];
1304:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
1305:       }
1306:     }
1307:     PetscFree2(xw,ww);
1308:     break;
1309:   case 3:
1310:     PetscMalloc2(npoints,&xw,npoints,&ww);
1311:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
1312:     for (i = 0; i < npoints; ++i) {
1313:       for (j = 0; j < npoints; ++j) {
1314:         for (k = 0; k < npoints; ++k) {
1315:           x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
1316:           x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
1317:           x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
1318:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
1319:         }
1320:       }
1321:     }
1322:     PetscFree2(xw,ww);
1323:     break;
1324:   default:
1325:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
1326:   }
1327:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
1328:   PetscQuadratureSetOrder(*q, 2*npoints-1);
1329:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
1330:   PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");
1331:   return(0);
1332: }

1334: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
1335: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
1336: {
1338:   *xi   = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
1339:   *eta  = 0.5  * (1.0 + y) * (1.0 - z) - 1.0;
1340:   *zeta = z;
1341:   return(0);
1342: }


1345: /*@
1346:   PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex

1348:   Not Collective

1350:   Input Arguments:
1351: + dim     - The simplex dimension
1352: . Nc      - The number of components
1353: . npoints - The number of points in one dimension
1354: . a       - left end of interval (often-1)
1355: - b       - right end of interval (often +1)

1357:   Output Argument:
1358: . q - A PetscQuadrature object

1360:   Level: intermediate

1362:   References:
1363: .  1. - Karniadakis and Sherwin.  FIAT

1365:   Note: For dim == 1, this is Gauss-Legendre quadrature

1367: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
1368: @*/
1369: PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1370: {
1371:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints;
1372:   PetscReal     *px, *wx, *py, *wy, *pz, *wz, *x, *w;

1376:   if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
1377:   PetscMalloc1(totpoints*dim, &x);
1378:   PetscMalloc1(totpoints*Nc, &w);
1379:   switch (dim) {
1380:   case 0:
1381:     PetscFree(x);
1382:     PetscFree(w);
1383:     PetscMalloc1(1, &x);
1384:     PetscMalloc1(Nc, &w);
1385:     x[0] = 0.0;
1386:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
1387:     break;
1388:   case 1:
1389:     PetscMalloc1(npoints,&wx);
1390:     PetscDTGaussJacobiQuadrature(npoints, -1., 1., 0.0, 0.0, x, wx);
1391:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i];
1392:     PetscFree(wx);
1393:     break;
1394:   case 2:
1395:     PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);
1396:     PetscDTGaussJacobiQuadrature(npoints, -1., 1., 0.0, 0.0, px, wx);
1397:     PetscDTGaussJacobiQuadrature(npoints, -1., 1., 1.0, 0.0, py, wy);
1398:     for (i = 0; i < npoints; ++i) {
1399:       for (j = 0; j < npoints; ++j) {
1400:         PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);
1401:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j];
1402:       }
1403:     }
1404:     PetscFree4(px,wx,py,wy);
1405:     break;
1406:   case 3:
1407:     PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);
1408:     PetscDTGaussJacobiQuadrature(npoints, -1., 1., 0.0, 0.0, px, wx);
1409:     PetscDTGaussJacobiQuadrature(npoints, -1., 1., 1.0, 0.0, py, wy);
1410:     PetscDTGaussJacobiQuadrature(npoints, -1., 1., 2.0, 0.0, pz, wz);
1411:     for (i = 0; i < npoints; ++i) {
1412:       for (j = 0; j < npoints; ++j) {
1413:         for (k = 0; k < npoints; ++k) {
1414:           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]);
1415:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k];
1416:         }
1417:       }
1418:     }
1419:     PetscFree6(px,wx,py,wy,pz,wz);
1420:     break;
1421:   default:
1422:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
1423:   }
1424:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
1425:   PetscQuadratureSetOrder(*q, 2*npoints-1);
1426:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
1427:   PetscObjectChangeTypeName((PetscObject)*q,"GaussJacobi");
1428:   return(0);
1429: }

1431: /*@
1432:   PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell

1434:   Not Collective

1436:   Input Arguments:
1437: + dim   - The cell dimension
1438: . level - The number of points in one dimension, 2^l
1439: . a     - left end of interval (often-1)
1440: - b     - right end of interval (often +1)

1442:   Output Argument:
1443: . q - A PetscQuadrature object

1445:   Level: intermediate

1447: .seealso: PetscDTGaussTensorQuadrature()
1448: @*/
1449: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
1450: {
1451:   const PetscInt  p     = 16;                        /* Digits of precision in the evaluation */
1452:   const PetscReal alpha = (b-a)/2.;                  /* Half-width of the integration interval */
1453:   const PetscReal beta  = (b+a)/2.;                  /* Center of the integration interval */
1454:   const PetscReal h     = PetscPowReal(2.0, -level); /* Step size, length between x_k */
1455:   PetscReal       xk;                                /* Quadrature point x_k on reference domain [-1, 1] */
1456:   PetscReal       wk    = 0.5*PETSC_PI;              /* Quadrature weight at x_k */
1457:   PetscReal      *x, *w;
1458:   PetscInt        K, k, npoints;
1459:   PetscErrorCode  ierr;

1462:   if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
1463:   if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
1464:   /* Find K such that the weights are < 32 digits of precision */
1465:   for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
1466:     wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
1467:   }
1468:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
1469:   PetscQuadratureSetOrder(*q, 2*K+1);
1470:   npoints = 2*K-1;
1471:   PetscMalloc1(npoints*dim, &x);
1472:   PetscMalloc1(npoints, &w);
1473:   /* Center term */
1474:   x[0] = beta;
1475:   w[0] = 0.5*alpha*PETSC_PI;
1476:   for (k = 1; k < K; ++k) {
1477:     wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1478:     xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
1479:     x[2*k-1] = -alpha*xk+beta;
1480:     w[2*k-1] = wk;
1481:     x[2*k+0] =  alpha*xk+beta;
1482:     w[2*k+0] = wk;
1483:   }
1484:   PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
1485:   return(0);
1486: }

1488: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1489: {
1490:   const PetscInt  p     = 16;        /* Digits of precision in the evaluation */
1491:   const PetscReal alpha = (b-a)/2.;  /* Half-width of the integration interval */
1492:   const PetscReal beta  = (b+a)/2.;  /* Center of the integration interval */
1493:   PetscReal       h     = 1.0;       /* Step size, length between x_k */
1494:   PetscInt        l     = 0;         /* Level of refinement, h = 2^{-l} */
1495:   PetscReal       osum  = 0.0;       /* Integral on last level */
1496:   PetscReal       psum  = 0.0;       /* Integral on the level before the last level */
1497:   PetscReal       sum;               /* Integral on current level */
1498:   PetscReal       yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1499:   PetscReal       lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1500:   PetscReal       wk;                /* Quadrature weight at x_k */
1501:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
1502:   PetscInt        d;                 /* Digits of precision in the integral */

1505:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1506:   /* Center term */
1507:   func(beta, &lval);
1508:   sum = 0.5*alpha*PETSC_PI*lval;
1509:   /* */
1510:   do {
1511:     PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
1512:     PetscInt  k = 1;

1514:     ++l;
1515:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1516:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1517:     psum = osum;
1518:     osum = sum;
1519:     h   *= 0.5;
1520:     sum *= 0.5;
1521:     do {
1522:       wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1523:       yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1524:       lx = -alpha*(1.0 - yk)+beta;
1525:       rx =  alpha*(1.0 - yk)+beta;
1526:       func(lx, &lval);
1527:       func(rx, &rval);
1528:       lterm   = alpha*wk*lval;
1529:       maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
1530:       sum    += lterm;
1531:       rterm   = alpha*wk*rval;
1532:       maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
1533:       sum    += rterm;
1534:       ++k;
1535:       /* Only need to evaluate every other point on refined levels */
1536:       if (l != 1) ++k;
1537:     } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */

1539:     d1 = PetscLog10Real(PetscAbsReal(sum - osum));
1540:     d2 = PetscLog10Real(PetscAbsReal(sum - psum));
1541:     d3 = PetscLog10Real(maxTerm) - p;
1542:     if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
1543:     else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
1544:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1545:   } while (d < digits && l < 12);
1546:   *sol = sum;

1548:   return(0);
1549: }

1551: #if defined(PETSC_HAVE_MPFR)
1552: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1553: {
1554:   const PetscInt  safetyFactor = 2;  /* Calculate abcissa until 2*p digits */
1555:   PetscInt        l            = 0;  /* Level of refinement, h = 2^{-l} */
1556:   mpfr_t          alpha;             /* Half-width of the integration interval */
1557:   mpfr_t          beta;              /* Center of the integration interval */
1558:   mpfr_t          h;                 /* Step size, length between x_k */
1559:   mpfr_t          osum;              /* Integral on last level */
1560:   mpfr_t          psum;              /* Integral on the level before the last level */
1561:   mpfr_t          sum;               /* Integral on current level */
1562:   mpfr_t          yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1563:   mpfr_t          lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1564:   mpfr_t          wk;                /* Quadrature weight at x_k */
1565:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
1566:   PetscInt        d;                 /* Digits of precision in the integral */
1567:   mpfr_t          pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;

1570:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1571:   /* Create high precision storage */
1572:   mpfr_inits2(PetscCeilReal(safetyFactor*digits*PetscLogReal(10.)/PetscLogReal(2.)), alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
1573:   /* Initialization */
1574:   mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
1575:   mpfr_set_d(beta,  0.5*(b+a), MPFR_RNDN);
1576:   mpfr_set_d(osum,  0.0,       MPFR_RNDN);
1577:   mpfr_set_d(psum,  0.0,       MPFR_RNDN);
1578:   mpfr_set_d(h,     1.0,       MPFR_RNDN);
1579:   mpfr_const_pi(pi2, MPFR_RNDN);
1580:   mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
1581:   /* Center term */
1582:   func(0.5*(b+a), &lval);
1583:   mpfr_set(sum, pi2, MPFR_RNDN);
1584:   mpfr_mul(sum, sum, alpha, MPFR_RNDN);
1585:   mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
1586:   /* */
1587:   do {
1588:     PetscReal d1, d2, d3, d4;
1589:     PetscInt  k = 1;

1591:     ++l;
1592:     mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
1593:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1594:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1595:     mpfr_set(psum, osum, MPFR_RNDN);
1596:     mpfr_set(osum,  sum, MPFR_RNDN);
1597:     mpfr_mul_d(h,   h,   0.5, MPFR_RNDN);
1598:     mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
1599:     do {
1600:       mpfr_set_si(kh, k, MPFR_RNDN);
1601:       mpfr_mul(kh, kh, h, MPFR_RNDN);
1602:       /* Weight */
1603:       mpfr_set(wk, h, MPFR_RNDN);
1604:       mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
1605:       mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
1606:       mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
1607:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
1608:       mpfr_sqr(tmp, tmp, MPFR_RNDN);
1609:       mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
1610:       mpfr_div(wk, wk, tmp, MPFR_RNDN);
1611:       /* Abscissa */
1612:       mpfr_set_d(yk, 1.0, MPFR_RNDZ);
1613:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
1614:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1615:       mpfr_exp(tmp, msinh, MPFR_RNDN);
1616:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1617:       /* Quadrature points */
1618:       mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
1619:       mpfr_mul(lx, lx, alpha, MPFR_RNDU);
1620:       mpfr_add(lx, lx, beta, MPFR_RNDU);
1621:       mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
1622:       mpfr_mul(rx, rx, alpha, MPFR_RNDD);
1623:       mpfr_add(rx, rx, beta, MPFR_RNDD);
1624:       /* Evaluation */
1625:       func(mpfr_get_d(lx, MPFR_RNDU), &lval);
1626:       func(mpfr_get_d(rx, MPFR_RNDD), &rval);
1627:       /* Update */
1628:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1629:       mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
1630:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
1631:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1632:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1633:       mpfr_set(curTerm, tmp, MPFR_RNDN);
1634:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1635:       mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
1636:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
1637:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1638:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1639:       mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
1640:       ++k;
1641:       /* Only need to evaluate every other point on refined levels */
1642:       if (l != 1) ++k;
1643:       mpfr_log10(tmp, wk, MPFR_RNDN);
1644:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1645:     } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
1646:     mpfr_sub(tmp, sum, osum, MPFR_RNDN);
1647:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1648:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1649:     d1 = mpfr_get_d(tmp, MPFR_RNDN);
1650:     mpfr_sub(tmp, sum, psum, MPFR_RNDN);
1651:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1652:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1653:     d2 = mpfr_get_d(tmp, MPFR_RNDN);
1654:     mpfr_log10(tmp, maxTerm, MPFR_RNDN);
1655:     d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
1656:     mpfr_log10(tmp, curTerm, MPFR_RNDN);
1657:     d4 = mpfr_get_d(tmp, MPFR_RNDN);
1658:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1659:   } while (d < digits && l < 8);
1660:   *sol = mpfr_get_d(sum, MPFR_RNDN);
1661:   /* Cleanup */
1662:   mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
1663:   return(0);
1664: }
1665: #else

1667: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1668: {
1669:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
1670: }
1671: #endif

1673: /* Overwrites A. Can only handle full-rank problems with m>=n
1674:  * A in column-major format
1675:  * Ainv in row-major format
1676:  * tau has length m
1677:  * worksize must be >= max(1,n)
1678:  */
1679: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1680: {
1682:   PetscBLASInt   M,N,K,lda,ldb,ldwork,info;
1683:   PetscScalar    *A,*Ainv,*R,*Q,Alpha;

1686: #if defined(PETSC_USE_COMPLEX)
1687:   {
1688:     PetscInt i,j;
1689:     PetscMalloc2(m*n,&A,m*n,&Ainv);
1690:     for (j=0; j<n; j++) {
1691:       for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1692:     }
1693:     mstride = m;
1694:   }
1695: #else
1696:   A = A_in;
1697:   Ainv = Ainv_out;
1698: #endif

1700:   PetscBLASIntCast(m,&M);
1701:   PetscBLASIntCast(n,&N);
1702:   PetscBLASIntCast(mstride,&lda);
1703:   PetscBLASIntCast(worksize,&ldwork);
1704:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1705:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1706:   PetscFPTrapPop();
1707:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1708:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1710:   /* Extract an explicit representation of Q */
1711:   Q = Ainv;
1712:   PetscArraycpy(Q,A,mstride*n);
1713:   K = N;                        /* full rank */
1714:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1715:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

1717:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1718:   Alpha = 1.0;
1719:   ldb = lda;
1720:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1721:   /* Ainv is Q, overwritten with inverse */

1723: #if defined(PETSC_USE_COMPLEX)
1724:   {
1725:     PetscInt i;
1726:     for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1727:     PetscFree2(A,Ainv);
1728:   }
1729: #endif
1730:   return(0);
1731: }

1733: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1734: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1735: {
1737:   PetscReal      *Bv;
1738:   PetscInt       i,j;

1741:   PetscMalloc1((ninterval+1)*ndegree,&Bv);
1742:   /* Point evaluation of L_p on all the source vertices */
1743:   PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1744:   /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1745:   for (i=0; i<ninterval; i++) {
1746:     for (j=0; j<ndegree; j++) {
1747:       if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1748:       else           B[i*ndegree+j]   = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1749:     }
1750:   }
1751:   PetscFree(Bv);
1752:   return(0);
1753: }

1755: /*@
1756:    PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals

1758:    Not Collective

1760:    Input Arguments:
1761: +  degree - degree of reconstruction polynomial
1762: .  nsource - number of source intervals
1763: .  sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1764: .  ntarget - number of target intervals
1765: -  targetx - sorted coordinates of target cell boundaries (length ntarget+1)

1767:    Output Arguments:
1768: .  R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]

1770:    Level: advanced

1772: .seealso: PetscDTLegendreEval()
1773: @*/
1774: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1775: {
1777:   PetscInt       i,j,k,*bdegrees,worksize;
1778:   PetscReal      xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1779:   PetscScalar    *tau,*work;

1785:   if (degree >= nsource) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Reconstruction degree %D must be less than number of source intervals %D",degree,nsource);
1786: #if defined(PETSC_USE_DEBUG)
1787:   for (i=0; i<nsource; i++) {
1788:     if (sourcex[i] >= sourcex[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Source interval %D has negative orientation (%g,%g)",i,(double)sourcex[i],(double)sourcex[i+1]);
1789:   }
1790:   for (i=0; i<ntarget; i++) {
1791:     if (targetx[i] >= targetx[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Target interval %D has negative orientation (%g,%g)",i,(double)targetx[i],(double)targetx[i+1]);
1792:   }
1793: #endif
1794:   xmin = PetscMin(sourcex[0],targetx[0]);
1795:   xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1796:   center = (xmin + xmax)/2;
1797:   hscale = (xmax - xmin)/2;
1798:   worksize = nsource;
1799:   PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1800:   PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1801:   for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1802:   for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1803:   PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1804:   PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1805:   for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1806:   PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1807:   for (i=0; i<ntarget; i++) {
1808:     PetscReal rowsum = 0;
1809:     for (j=0; j<nsource; j++) {
1810:       PetscReal sum = 0;
1811:       for (k=0; k<degree+1; k++) {
1812:         sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1813:       }
1814:       R[i*nsource+j] = sum;
1815:       rowsum += sum;
1816:     }
1817:     for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1818:   }
1819:   PetscFree4(bdegrees,sourcey,Bsource,work);
1820:   PetscFree4(tau,Bsinv,targety,Btarget);
1821:   return(0);
1822: }

1824: /*@C
1825:    PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points

1827:    Not Collective

1829:    Input Parameter:
1830: +  n - the number of GLL nodes
1831: .  nodes - the GLL nodes
1832: .  weights - the GLL weights
1833: -  f - the function values at the nodes

1835:    Output Parameter:
1836: .  in - the value of the integral

1838:    Level: beginner

1840: .seealso: PetscDTGaussLobattoLegendreQuadrature()

1842: @*/
1843: PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal *nodes,PetscReal *weights,const PetscReal *f,PetscReal *in)
1844: {
1845:   PetscInt          i;

1848:   *in = 0.;
1849:   for (i=0; i<n; i++) {
1850:     *in += f[i]*f[i]*weights[i];
1851:   }
1852:   return(0);
1853: }

1855: /*@C
1856:    PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element

1858:    Not Collective

1860:    Input Parameter:
1861: +  n - the number of GLL nodes
1862: .  nodes - the GLL nodes
1863: -  weights - the GLL weights

1865:    Output Parameter:
1866: .  A - the stiffness element

1868:    Level: beginner

1870:    Notes:
1871:     Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy()

1873:    You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented (the array is symmetric)

1875: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()

1877: @*/
1878: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1879: {
1880:   PetscReal        **A;
1881:   PetscErrorCode  ierr;
1882:   const PetscReal  *gllnodes = nodes;
1883:   const PetscInt   p = n-1;
1884:   PetscReal        z0,z1,z2 = -1,x,Lpj,Lpr;
1885:   PetscInt         i,j,nn,r;

1888:   PetscMalloc1(n,&A);
1889:   PetscMalloc1(n*n,&A[0]);
1890:   for (i=1; i<n; i++) A[i] = A[i-1]+n;

1892:   for (j=1; j<p; j++) {
1893:     x  = gllnodes[j];
1894:     z0 = 1.;
1895:     z1 = x;
1896:     for (nn=1; nn<p; nn++) {
1897:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1898:       z0 = z1;
1899:       z1 = z2;
1900:     }
1901:     Lpj=z2;
1902:     for (r=1; r<p; r++) {
1903:       if (r == j) {
1904:         A[j][j]=2./(3.*(1.-gllnodes[j]*gllnodes[j])*Lpj*Lpj);
1905:       } else {
1906:         x  = gllnodes[r];
1907:         z0 = 1.;
1908:         z1 = x;
1909:         for (nn=1; nn<p; nn++) {
1910:           z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1911:           z0 = z1;
1912:           z1 = z2;
1913:         }
1914:         Lpr     = z2;
1915:         A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(gllnodes[j]-gllnodes[r])*(gllnodes[j]-gllnodes[r]));
1916:       }
1917:     }
1918:   }
1919:   for (j=1; j<p+1; j++) {
1920:     x  = gllnodes[j];
1921:     z0 = 1.;
1922:     z1 = x;
1923:     for (nn=1; nn<p; nn++) {
1924:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1925:       z0 = z1;
1926:       z1 = z2;
1927:     }
1928:     Lpj     = z2;
1929:     A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+gllnodes[j])*(1.+gllnodes[j]));
1930:     A[0][j] = A[j][0];
1931:   }
1932:   for (j=0; j<p; j++) {
1933:     x  = gllnodes[j];
1934:     z0 = 1.;
1935:     z1 = x;
1936:     for (nn=1; nn<p; nn++) {
1937:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1938:       z0 = z1;
1939:       z1 = z2;
1940:     }
1941:     Lpj=z2;

1943:     A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-gllnodes[j])*(1.-gllnodes[j]));
1944:     A[j][p] = A[p][j];
1945:   }
1946:   A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.;
1947:   A[p][p]=A[0][0];
1948:   *AA = A;
1949:   return(0);
1950: }

1952: /*@C
1953:    PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element

1955:    Not Collective

1957:    Input Parameter:
1958: +  n - the number of GLL nodes
1959: .  nodes - the GLL nodes
1960: .  weights - the GLL weightss
1961: -  A - the stiffness element

1963:    Level: beginner

1965: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate()

1967: @*/
1968: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1969: {

1973:   PetscFree((*AA)[0]);
1974:   PetscFree(*AA);
1975:   *AA  = NULL;
1976:   return(0);
1977: }

1979: /*@C
1980:    PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element

1982:    Not Collective

1984:    Input Parameter:
1985: +  n - the number of GLL nodes
1986: .  nodes - the GLL nodes
1987: .  weights - the GLL weights

1989:    Output Parameter:
1990: .  AA - the stiffness element
1991: -  AAT - the transpose of AA (pass in NULL if you do not need this array)

1993:    Level: beginner

1995:    Notes:
1996:     Destroy this with PetscGaussLobattoLegendreElementGradientDestroy()

1998:    You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented

2000: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()

2002: @*/
2003: PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
2004: {
2005:   PetscReal        **A, **AT = NULL;
2006:   PetscErrorCode  ierr;
2007:   const PetscReal  *gllnodes = nodes;
2008:   const PetscInt   p = n-1;
2009:   PetscReal        Li, Lj,d0;
2010:   PetscInt         i,j;

2013:   PetscMalloc1(n,&A);
2014:   PetscMalloc1(n*n,&A[0]);
2015:   for (i=1; i<n; i++) A[i] = A[i-1]+n;

2017:   if (AAT) {
2018:     PetscMalloc1(n,&AT);
2019:     PetscMalloc1(n*n,&AT[0]);
2020:     for (i=1; i<n; i++) AT[i] = AT[i-1]+n;
2021:   }

2023:   if (n==1) {A[0][0] = 0.;}
2024:   d0 = (PetscReal)p*((PetscReal)p+1.)/4.;
2025:   for  (i=0; i<n; i++) {
2026:     for  (j=0; j<n; j++) {
2027:       A[i][j] = 0.;
2028:       PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li);
2029:       PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj);
2030:       if (i!=j)             A[i][j] = Li/(Lj*(gllnodes[i]-gllnodes[j]));
2031:       if ((j==i) && (i==0)) A[i][j] = -d0;
2032:       if (j==i && i==p)     A[i][j] = d0;
2033:       if (AT) AT[j][i] = A[i][j];
2034:     }
2035:   }
2036:   if (AAT) *AAT = AT;
2037:   *AA  = A;
2038:   return(0);
2039: }

2041: /*@C
2042:    PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGaussLobattoLegendreElementGradientCreate()

2044:    Not Collective

2046:    Input Parameter:
2047: +  n - the number of GLL nodes
2048: .  nodes - the GLL nodes
2049: .  weights - the GLL weights
2050: .  AA - the stiffness element
2051: -  AAT - the transpose of the element

2053:    Level: beginner

2055: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionCreate()

2057: @*/
2058: PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
2059: {

2063:   PetscFree((*AA)[0]);
2064:   PetscFree(*AA);
2065:   *AA  = NULL;
2066:   if (*AAT) {
2067:     PetscFree((*AAT)[0]);
2068:     PetscFree(*AAT);
2069:     *AAT  = NULL;
2070:   }
2071:   return(0);
2072: }

2074: /*@C
2075:    PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element

2077:    Not Collective

2079:    Input Parameter:
2080: +  n - the number of GLL nodes
2081: .  nodes - the GLL nodes
2082: -  weights - the GLL weightss

2084:    Output Parameter:
2085: .  AA - the stiffness element

2087:    Level: beginner

2089:    Notes:
2090:     Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy()

2092:    This is the same as the Gradient operator multiplied by the diagonal mass matrix

2094:    You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented

2096: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionDestroy()

2098: @*/
2099: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2100: {
2101:   PetscReal       **D;
2102:   PetscErrorCode  ierr;
2103:   const PetscReal  *gllweights = weights;
2104:   const PetscInt   glln = n;
2105:   PetscInt         i,j;

2108:   PetscGaussLobattoLegendreElementGradientCreate(n,nodes,weights,&D,NULL);
2109:   for (i=0; i<glln; i++){
2110:     for (j=0; j<glln; j++) {
2111:       D[i][j] = gllweights[i]*D[i][j];
2112:     }
2113:   }
2114:   *AA = D;
2115:   return(0);
2116: }

2118: /*@C
2119:    PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element

2121:    Not Collective

2123:    Input Parameter:
2124: +  n - the number of GLL nodes
2125: .  nodes - the GLL nodes
2126: .  weights - the GLL weights
2127: -  A - advection

2129:    Level: beginner

2131: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementAdvectionCreate()

2133: @*/
2134: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2135: {

2139:   PetscFree((*AA)[0]);
2140:   PetscFree(*AA);
2141:   *AA  = NULL;
2142:   return(0);
2143: }

2145: PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2146: {
2147:   PetscReal        **A;
2148:   PetscErrorCode  ierr;
2149:   const PetscReal  *gllweights = weights;
2150:   const PetscInt   glln = n;
2151:   PetscInt         i,j;

2154:   PetscMalloc1(glln,&A);
2155:   PetscMalloc1(glln*glln,&A[0]);
2156:   for (i=1; i<glln; i++) A[i] = A[i-1]+glln;
2157:   if (glln==1) {A[0][0] = 0.;}
2158:   for  (i=0; i<glln; i++) {
2159:     for  (j=0; j<glln; j++) {
2160:       A[i][j] = 0.;
2161:       if (j==i)     A[i][j] = gllweights[i];
2162:     }
2163:   }
2164:   *AA  = A;
2165:   return(0);
2166: }

2168: PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2169: {

2173:   PetscFree((*AA)[0]);
2174:   PetscFree(*AA);
2175:   *AA  = NULL;
2176:   return(0);
2177: }

2179: /*@
2180:   PetscDTIndexToBary - convert an index into a barycentric coordinate.

2182:   Input Parameters:
2183: + len - the desired length of the barycentric tuple (usually 1 more than the dimension it represents, so a barycentric coordinate in a triangle has length 3)
2184: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
2185: - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum)

2187:   Output Parameter:
2188: . coord - will be filled with the barycentric coordinate

2190:   Level: beginner

2192:   Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the
2193:   least significant and the last index is the most significant.

2195: .seealso: PetscDTBaryToIndex
2196: @*/
2197: PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[])
2198: {
2199:   PetscInt c, d, s, total, subtotal, nexttotal;

2202:   if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
2203:   if (index < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative");
2204:   if (!len) {
2205:     if (!sum && !index) return(0);
2206:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
2207:   }
2208:   for (c = 1, total = 1; c <= len; c++) {
2209:     /* total is the number of ways to have a tuple of length c with sum */
2210:     if (index < total) break;
2211:     total = (total * (sum + c)) / c;
2212:   }
2213:   if (c > len) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range");
2214:   for (d = c; d < len; d++) coord[d] = 0;
2215:   for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) {
2216:     /* subtotal is the number of ways to have a tuple of length c with sum s */
2217:     /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */
2218:     if ((index + subtotal) >= total) {
2219:       coord[--c] = sum - s;
2220:       index -= (total - subtotal);
2221:       sum = s;
2222:       total = nexttotal;
2223:       subtotal = 1;
2224:       nexttotal = 1;
2225:       s = 0;
2226:     } else {
2227:       subtotal = (subtotal * (c + s)) / (s + 1);
2228:       nexttotal = (nexttotal * (c - 1 + s)) / (s + 1);
2229:       s++;
2230:     }
2231:   }
2232:   return(0);
2233: }

2235: /*@
2236:   PetscDTBaryToIndex - convert a barycentric coordinate to an index

2238:   Input Parameters:
2239: + len - the desired length of the barycentric tuple (usually 1 more than the dimension it represents, so a barycentric coordinate in a triangle has length 3)
2240: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
2241: - coord - a barycentric coordinate with the given length and sum

2243:   Output Parameter:
2244: . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum)

2246:   Level: beginner

2248:   Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the
2249:   least significant and the last index is the most significant.

2251: .seealso: PetscDTIndexToBary
2252: @*/
2253: PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index)
2254: {
2255:   PetscInt c;
2256:   PetscInt i;
2257:   PetscInt total;

2260:   if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
2261:   if (!len) {
2262:     if (!sum) {
2263:       *index = 0;
2264:       return(0);
2265:     }
2266:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
2267:   }
2268:   for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c;
2269:   i = total - 1;
2270:   c = len - 1;
2271:   sum -= coord[c];
2272:   while (sum > 0) {
2273:     PetscInt subtotal;
2274:     PetscInt s;

2276:     for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s;
2277:     i   -= subtotal;
2278:     sum -= coord[--c];
2279:   }
2280:   *index = i;
2281:   return(0);
2282: }