Actual source code: dt.c

petsc-3.11.1 2019-04-17
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: 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";

 25: /*@
 26:   PetscQuadratureCreate - Create a PetscQuadrature object

 28:   Collective on MPI_Comm

 30:   Input Parameter:
 31: . comm - The communicator for the PetscQuadrature object

 33:   Output Parameter:
 34: . q  - The PetscQuadrature object

 36:   Level: beginner

 38: .keywords: PetscQuadrature, quadrature, create
 39: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
 40: @*/
 41: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
 42: {

 47:   PetscSysInitializePackage();
 48:   PetscHeaderCreate(*q,PETSC_OBJECT_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
 49:   (*q)->dim       = -1;
 50:   (*q)->Nc        =  1;
 51:   (*q)->order     = -1;
 52:   (*q)->numPoints = 0;
 53:   (*q)->points    = NULL;
 54:   (*q)->weights   = NULL;
 55:   return(0);
 56: }

 58: /*@
 59:   PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object

 61:   Collective on PetscQuadrature

 63:   Input Parameter:
 64: . q  - The PetscQuadrature object

 66:   Output Parameter:
 67: . r  - The new PetscQuadrature object

 69:   Level: beginner

 71: .keywords: PetscQuadrature, quadrature, clone
 72: .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData()
 73: @*/
 74: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
 75: {
 76:   PetscInt         order, dim, Nc, Nq;
 77:   const PetscReal *points, *weights;
 78:   PetscReal       *p, *w;
 79:   PetscErrorCode   ierr;

 83:   PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);
 84:   PetscQuadratureGetOrder(q, &order);
 85:   PetscQuadratureSetOrder(*r, order);
 86:   PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);
 87:   PetscMalloc1(Nq*dim, &p);
 88:   PetscMalloc1(Nq*Nc, &w);
 89:   PetscMemcpy(p, points, Nq*dim * sizeof(PetscReal));
 90:   PetscMemcpy(w, weights, Nc * Nq * sizeof(PetscReal));
 91:   PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);
 92:   return(0);
 93: }

 95: /*@
 96:   PetscQuadratureDestroy - Destroys a PetscQuadrature object

 98:   Collective on PetscQuadrature

100:   Input Parameter:
101: . q  - The PetscQuadrature object

103:   Level: beginner

105: .keywords: PetscQuadrature, quadrature, destroy
106: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
107: @*/
108: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
109: {

113:   if (!*q) return(0);
115:   if (--((PetscObject)(*q))->refct > 0) {
116:     *q = NULL;
117:     return(0);
118:   }
119:   PetscFree((*q)->points);
120:   PetscFree((*q)->weights);
121:   PetscHeaderDestroy(q);
122:   return(0);
123: }

125: /*@
126:   PetscQuadratureGetOrder - Return the order of the method

128:   Not collective

130:   Input Parameter:
131: . q - The PetscQuadrature object

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

136:   Level: intermediate

138: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
139: @*/
140: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
141: {
145:   *order = q->order;
146:   return(0);
147: }

149: /*@
150:   PetscQuadratureSetOrder - Return the order of the method

152:   Not collective

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

158:   Level: intermediate

160: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
161: @*/
162: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
163: {
166:   q->order = order;
167:   return(0);
168: }

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

173:   Not collective

175:   Input Parameter:
176: . q - The PetscQuadrature object

178:   Output Parameter:
179: . Nc - The number of components

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

183:   Level: intermediate

185: .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
186: @*/
187: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
188: {
192:   *Nc = q->Nc;
193:   return(0);
194: }

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

199:   Not collective

201:   Input Parameters:
202: + q  - The PetscQuadrature object
203: - Nc - The number of components

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

207:   Level: intermediate

209: .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
210: @*/
211: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
212: {
215:   q->Nc = Nc;
216:   return(0);
217: }

219: /*@C
220:   PetscQuadratureGetData - Returns the data defining the quadrature

222:   Not collective

224:   Input Parameter:
225: . q  - The PetscQuadrature object

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

234:   Level: intermediate

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

239: .keywords: PetscQuadrature, quadrature
240: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
241: @*/
242: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
243: {
246:   if (dim) {
248:     *dim = q->dim;
249:   }
250:   if (Nc) {
252:     *Nc = q->Nc;
253:   }
254:   if (npoints) {
256:     *npoints = q->numPoints;
257:   }
258:   if (points) {
260:     *points = q->points;
261:   }
262:   if (weights) {
264:     *weights = q->weights;
265:   }
266:   return(0);
267: }

269: /*@C
270:   PetscQuadratureSetData - Sets the data defining the quadrature

272:   Not collective

274:   Input Parameters:
275: + q  - The PetscQuadrature object
276: . dim - The spatial dimension
277: . Nc - The number of components
278: . npoints - The number of quadrature points
279: . points - The coordinates of each quadrature point
280: - weights - The weight of each quadrature point

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

284:   Level: intermediate

286: .keywords: PetscQuadrature, quadrature
287: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
288: @*/
289: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
290: {
293:   if (dim >= 0)     q->dim       = dim;
294:   if (Nc >= 0)      q->Nc        = Nc;
295:   if (npoints >= 0) q->numPoints = npoints;
296:   if (points) {
298:     q->points = points;
299:   }
300:   if (weights) {
302:     q->weights = weights;
303:   }
304:   return(0);
305: }

307: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
308: {
309:   PetscInt          q, d, c;
310:   PetscViewerFormat format;
311:   PetscErrorCode    ierr;

314:   if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "Quadrature on %D points with %D components of order %D\n", quad->numPoints, quad->Nc, quad->order);}
315:   else              {PetscViewerASCIIPrintf(v, "Quadrature on %D points of order %D\n", quad->numPoints, quad->order);}
316:   PetscViewerGetFormat(v, &format);
317:   if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
318:   for (q = 0; q < quad->numPoints; ++q) {
319:     PetscViewerASCIIPrintf(v, "(");
320:     PetscViewerASCIIUseTabs(v, PETSC_FALSE);
321:     for (d = 0; d < quad->dim; ++d) {
322:       if (d) {PetscViewerASCIIPrintf(v, ", ");}
323:       PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q*quad->dim+d]);
324:     }
325:     PetscViewerASCIIPrintf(v, ") ");
326:     if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "(");}
327:     for (c = 0; c < quad->Nc; ++c) {
328:       if (c) {PetscViewerASCIIPrintf(v, ", ");}
329:       PetscViewerASCIIPrintf(v, "%g", (double)quad->weights[q*quad->Nc+c]);
330:     }
331:     if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, ")");}
332:     PetscViewerASCIIPrintf(v, "\n");
333:     PetscViewerASCIIUseTabs(v, PETSC_TRUE);
334:   }
335:   return(0);
336: }

338: /*@C
339:   PetscQuadratureView - Views a PetscQuadrature object

341:   Collective on PetscQuadrature

343:   Input Parameters:
344: + quad  - The PetscQuadrature object
345: - viewer - The PetscViewer object

347:   Level: beginner

349: .keywords: PetscQuadrature, quadrature, view
350: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
351: @*/
352: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
353: {
354:   PetscBool      iascii;

360:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) quad), &viewer);}
361:   PetscObjectPrintClassNamePrefixType((PetscObject)quad, viewer);
362:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
363:   PetscViewerASCIIPushTab(viewer);
364:   if (iascii) {PetscQuadratureView_Ascii(quad, viewer);}
365:   PetscViewerASCIIPopTab(viewer);
366:   return(0);
367: }

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

372:   Not collective

374:   Input Parameter:
375: + q - The original PetscQuadrature
376: . numSubelements - The number of subelements the original element is divided into
377: . v0 - An array of the initial points for each subelement
378: - jac - An array of the Jacobian mappings from the reference to each subelement

380:   Output Parameters:
381: . dim - The dimension

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

385:  Not available from Fortran

387:   Level: intermediate

389: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
390: @*/
391: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
392: {
393:   const PetscReal *points,    *weights;
394:   PetscReal       *pointsRef, *weightsRef;
395:   PetscInt         dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
396:   PetscErrorCode   ierr;

403:   PetscQuadratureCreate(PETSC_COMM_SELF, qref);
404:   PetscQuadratureGetOrder(q, &order);
405:   PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
406:   npointsRef = npoints*numSubelements;
407:   PetscMalloc1(npointsRef*dim,&pointsRef);
408:   PetscMalloc1(npointsRef*Nc, &weightsRef);
409:   for (c = 0; c < numSubelements; ++c) {
410:     for (p = 0; p < npoints; ++p) {
411:       for (d = 0; d < dim; ++d) {
412:         pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
413:         for (e = 0; e < dim; ++e) {
414:           pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
415:         }
416:       }
417:       /* Could also use detJ here */
418:       for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
419:     }
420:   }
421:   PetscQuadratureSetOrder(*qref, order);
422:   PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
423:   return(0);
424: }

426: /*@
427:    PetscDTLegendreEval - evaluate Legendre polynomial at points

429:    Not Collective

431:    Input Arguments:
432: +  npoints - number of spatial points to evaluate at
433: .  points - array of locations to evaluate at
434: .  ndegree - number of basis degrees to evaluate
435: -  degrees - sorted array of degrees to evaluate

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

442:    Level: intermediate

444: .seealso: PetscDTGaussQuadrature()
445: @*/
446: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
447: {
448:   PetscInt i,maxdegree;

451:   if (!npoints || !ndegree) return(0);
452:   maxdegree = degrees[ndegree-1];
453:   for (i=0; i<npoints; i++) {
454:     PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
455:     PetscInt  j,k;
456:     x    = points[i];
457:     pm2  = 0;
458:     pm1  = 1;
459:     pd2  = 0;
460:     pd1  = 0;
461:     pdd2 = 0;
462:     pdd1 = 0;
463:     k    = 0;
464:     if (degrees[k] == 0) {
465:       if (B) B[i*ndegree+k] = pm1;
466:       if (D) D[i*ndegree+k] = pd1;
467:       if (D2) D2[i*ndegree+k] = pdd1;
468:       k++;
469:     }
470:     for (j=1; j<=maxdegree; j++,k++) {
471:       PetscReal p,d,dd;
472:       p    = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
473:       d    = pd2 + (2*j-1)*pm1;
474:       dd   = pdd2 + (2*j-1)*pd1;
475:       pm2  = pm1;
476:       pm1  = p;
477:       pd2  = pd1;
478:       pd1  = d;
479:       pdd2 = pdd1;
480:       pdd1 = dd;
481:       if (degrees[k] == j) {
482:         if (B) B[i*ndegree+k] = p;
483:         if (D) D[i*ndegree+k] = d;
484:         if (D2) D2[i*ndegree+k] = dd;
485:       }
486:     }
487:   }
488:   return(0);
489: }

491: /*@
492:    PetscDTGaussQuadrature - create Gauss quadrature

494:    Not Collective

496:    Input Arguments:
497: +  npoints - number of points
498: .  a - left end of interval (often-1)
499: -  b - right end of interval (often +1)

501:    Output Arguments:
502: +  x - quadrature points
503: -  w - quadrature weights

505:    Level: intermediate

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

510: .seealso: PetscDTLegendreEval()
511: @*/
512: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
513: {
515:   PetscInt       i;
516:   PetscReal      *work;
517:   PetscScalar    *Z;
518:   PetscBLASInt   N,LDZ,info;

521:   PetscCitationsRegister(GaussCitation, &GaussCite);
522:   /* Set up the Golub-Welsch system */
523:   for (i=0; i<npoints; i++) {
524:     x[i] = 0;                   /* diagonal is 0 */
525:     if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
526:   }
527:   PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
528:   PetscBLASIntCast(npoints,&N);
529:   LDZ  = N;
530:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
531:   PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
532:   PetscFPTrapPop();
533:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");

535:   for (i=0; i<(npoints+1)/2; i++) {
536:     PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
537:     x[i]           = (a+b)/2 - y*(b-a)/2;
538:     if (x[i] == -0.0) x[i] = 0.0;
539:     x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;

541:     w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
542:   }
543:   PetscFree2(Z,work);
544:   return(0);
545: }

547: /*@
548:   PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature

550:   Not Collective

552:   Input Arguments:
553: + dim     - The spatial dimension
554: . Nc      - The number of components
555: . npoints - number of points in one dimension
556: . a       - left end of interval (often-1)
557: - b       - right end of interval (often +1)

559:   Output Argument:
560: . q - A PetscQuadrature object

562:   Level: intermediate

564: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
565: @*/
566: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
567: {
568:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
569:   PetscReal     *x, *w, *xw, *ww;

573:   PetscMalloc1(totpoints*dim,&x);
574:   PetscMalloc1(totpoints*Nc,&w);
575:   /* Set up the Golub-Welsch system */
576:   switch (dim) {
577:   case 0:
578:     PetscFree(x);
579:     PetscFree(w);
580:     PetscMalloc1(1, &x);
581:     PetscMalloc1(Nc, &w);
582:     x[0] = 0.0;
583:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
584:     break;
585:   case 1:
586:     PetscMalloc1(npoints,&ww);
587:     PetscDTGaussQuadrature(npoints, a, b, x, ww);
588:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
589:     PetscFree(ww);
590:     break;
591:   case 2:
592:     PetscMalloc2(npoints,&xw,npoints,&ww);
593:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
594:     for (i = 0; i < npoints; ++i) {
595:       for (j = 0; j < npoints; ++j) {
596:         x[(i*npoints+j)*dim+0] = xw[i];
597:         x[(i*npoints+j)*dim+1] = xw[j];
598:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
599:       }
600:     }
601:     PetscFree2(xw,ww);
602:     break;
603:   case 3:
604:     PetscMalloc2(npoints,&xw,npoints,&ww);
605:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
606:     for (i = 0; i < npoints; ++i) {
607:       for (j = 0; j < npoints; ++j) {
608:         for (k = 0; k < npoints; ++k) {
609:           x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
610:           x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
611:           x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
612:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
613:         }
614:       }
615:     }
616:     PetscFree2(xw,ww);
617:     break;
618:   default:
619:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
620:   }
621:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
622:   PetscQuadratureSetOrder(*q, 2*npoints-1);
623:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
624:   PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");
625:   return(0);
626: }

628: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
629:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
630: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
631: {
632:   PetscReal f = 1.0;
633:   PetscInt  i;

636:   for (i = 1; i < n+1; ++i) f *= i;
637:   *factorial = f;
638:   return(0);
639: }

641: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
642:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
643: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
644: {
645:   PetscReal apb, pn1, pn2;
646:   PetscInt  k;

649:   if (!n) {*P = 1.0; return(0);}
650:   if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
651:   apb = a + b;
652:   pn2 = 1.0;
653:   pn1 = 0.5 * (a - b + (apb + 2.0) * x);
654:   *P  = 0.0;
655:   for (k = 2; k < n+1; ++k) {
656:     PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
657:     PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
658:     PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
659:     PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);

661:     a2  = a2 / a1;
662:     a3  = a3 / a1;
663:     a4  = a4 / a1;
664:     *P  = (a2 + a3 * x) * pn1 - a4 * pn2;
665:     pn2 = pn1;
666:     pn1 = *P;
667:   }
668:   return(0);
669: }

671: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
672: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
673: {
674:   PetscReal      nP;

678:   if (!n) {*P = 0.0; return(0);}
679:   PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
680:   *P   = 0.5 * (a + b + n + 1) * nP;
681:   return(0);
682: }

684: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
685: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
686: {
688:   *xi  = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
689:   *eta = y;
690:   return(0);
691: }

693: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
694: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
695: {
697:   *xi   = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
698:   *eta  = 0.5  * (1.0 + y) * (1.0 - z) - 1.0;
699:   *zeta = z;
700:   return(0);
701: }

703: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
704: {
705:   PetscInt       maxIter = 100;
706:   PetscReal      eps     = 1.0e-8;
707:   PetscReal      a1, a2, a3, a4, a5, a6;
708:   PetscInt       k;


713:   a1      = PetscPowReal(2.0, a+b+1);
714: #if defined(PETSC_HAVE_TGAMMA)
715:   a2      = PetscTGamma(a + npoints + 1);
716:   a3      = PetscTGamma(b + npoints + 1);
717:   a4      = PetscTGamma(a + b + npoints + 1);
718: #else
719:   {
720:     PetscInt ia, ib;

722:     ia = (PetscInt) a;
723:     ib = (PetscInt) b;
724:     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 */
725:       PetscDTFactorial_Internal(ia + npoints, &a2);
726:       PetscDTFactorial_Internal(ib + npoints, &a3);
727:       PetscDTFactorial_Internal(ia + ib + npoints, &a4);
728:     } else {
729:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
730:     }
731:   }
732: #endif

734:   PetscDTFactorial_Internal(npoints, &a5);
735:   a6   = a1 * a2 * a3 / a4 / a5;
736:   /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
737:    Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
738:   for (k = 0; k < npoints; ++k) {
739:     PetscReal r = -PetscCosReal((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP;
740:     PetscInt  j;

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

747:       for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
748:       PetscDTComputeJacobi(a, b, npoints, r, &f);
749:       PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
750:       delta = f / (fp - f * s);
751:       r     = r - delta;
752:       if (PetscAbsReal(delta) < eps) break;
753:     }
754:     x[k] = r;
755:     PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
756:     w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
757:   }
758:   return(0);
759: }

761: /*@
762:   PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex

764:   Not Collective

766:   Input Arguments:
767: + dim     - The simplex dimension
768: . Nc      - The number of components
769: . npoints - The number of points in one dimension
770: . a       - left end of interval (often-1)
771: - b       - right end of interval (often +1)

773:   Output Argument:
774: . q - A PetscQuadrature object

776:   Level: intermediate

778:   References:
779: .  1. - Karniadakis and Sherwin.  FIAT

781: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
782: @*/
783: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
784: {
785:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints;
786:   PetscReal     *px, *wx, *py, *wy, *pz, *wz, *x, *w;
787:   PetscInt       i, j, k, c;

791:   if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
792:   PetscMalloc1(totpoints*dim, &x);
793:   PetscMalloc1(totpoints*Nc, &w);
794:   switch (dim) {
795:   case 0:
796:     PetscFree(x);
797:     PetscFree(w);
798:     PetscMalloc1(1, &x);
799:     PetscMalloc1(Nc, &w);
800:     x[0] = 0.0;
801:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
802:     break;
803:   case 1:
804:     PetscMalloc1(npoints,&wx);
805:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, wx);
806:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i];
807:     PetscFree(wx);
808:     break;
809:   case 2:
810:     PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);
811:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
812:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
813:     for (i = 0; i < npoints; ++i) {
814:       for (j = 0; j < npoints; ++j) {
815:         PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);
816:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j];
817:       }
818:     }
819:     PetscFree4(px,wx,py,wy);
820:     break;
821:   case 3:
822:     PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);
823:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
824:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
825:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);
826:     for (i = 0; i < npoints; ++i) {
827:       for (j = 0; j < npoints; ++j) {
828:         for (k = 0; k < npoints; ++k) {
829:           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]);
830:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k];
831:         }
832:       }
833:     }
834:     PetscFree6(px,wx,py,wy,pz,wz);
835:     break;
836:   default:
837:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
838:   }
839:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
840:   PetscQuadratureSetOrder(*q, 2*npoints-1);
841:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
842:   PetscObjectChangeTypeName((PetscObject)*q,"GaussJacobi");
843:   return(0);
844: }

846: /*@
847:   PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell

849:   Not Collective

851:   Input Arguments:
852: + dim   - The cell dimension
853: . level - The number of points in one dimension, 2^l
854: . a     - left end of interval (often-1)
855: - b     - right end of interval (often +1)

857:   Output Argument:
858: . q - A PetscQuadrature object

860:   Level: intermediate

862: .seealso: PetscDTGaussTensorQuadrature()
863: @*/
864: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
865: {
866:   const PetscInt  p     = 16;                        /* Digits of precision in the evaluation */
867:   const PetscReal alpha = (b-a)/2.;                  /* Half-width of the integration interval */
868:   const PetscReal beta  = (b+a)/2.;                  /* Center of the integration interval */
869:   const PetscReal h     = PetscPowReal(2.0, -level); /* Step size, length between x_k */
870:   PetscReal       xk;                                /* Quadrature point x_k on reference domain [-1, 1] */
871:   PetscReal       wk    = 0.5*PETSC_PI;              /* Quadrature weight at x_k */
872:   PetscReal      *x, *w;
873:   PetscInt        K, k, npoints;
874:   PetscErrorCode  ierr;

877:   if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
878:   if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
879:   /* Find K such that the weights are < 32 digits of precision */
880:   for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
881:     wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
882:   }
883:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
884:   PetscQuadratureSetOrder(*q, 2*K+1);
885:   npoints = 2*K-1;
886:   PetscMalloc1(npoints*dim, &x);
887:   PetscMalloc1(npoints, &w);
888:   /* Center term */
889:   x[0] = beta;
890:   w[0] = 0.5*alpha*PETSC_PI;
891:   for (k = 1; k < K; ++k) {
892:     wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
893:     xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
894:     x[2*k-1] = -alpha*xk+beta;
895:     w[2*k-1] = wk;
896:     x[2*k+0] =  alpha*xk+beta;
897:     w[2*k+0] = wk;
898:   }
899:   PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
900:   return(0);
901: }

903: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
904: {
905:   const PetscInt  p     = 16;        /* Digits of precision in the evaluation */
906:   const PetscReal alpha = (b-a)/2.;  /* Half-width of the integration interval */
907:   const PetscReal beta  = (b+a)/2.;  /* Center of the integration interval */
908:   PetscReal       h     = 1.0;       /* Step size, length between x_k */
909:   PetscInt        l     = 0;         /* Level of refinement, h = 2^{-l} */
910:   PetscReal       osum  = 0.0;       /* Integral on last level */
911:   PetscReal       psum  = 0.0;       /* Integral on the level before the last level */
912:   PetscReal       sum;               /* Integral on current level */
913:   PetscReal       yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
914:   PetscReal       lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
915:   PetscReal       wk;                /* Quadrature weight at x_k */
916:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
917:   PetscInt        d;                 /* Digits of precision in the integral */

920:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
921:   /* Center term */
922:   func(beta, &lval);
923:   sum = 0.5*alpha*PETSC_PI*lval;
924:   /* */
925:   do {
926:     PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
927:     PetscInt  k = 1;

929:     ++l;
930:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
931:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
932:     psum = osum;
933:     osum = sum;
934:     h   *= 0.5;
935:     sum *= 0.5;
936:     do {
937:       wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
938:       yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
939:       lx = -alpha*(1.0 - yk)+beta;
940:       rx =  alpha*(1.0 - yk)+beta;
941:       func(lx, &lval);
942:       func(rx, &rval);
943:       lterm   = alpha*wk*lval;
944:       maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
945:       sum    += lterm;
946:       rterm   = alpha*wk*rval;
947:       maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
948:       sum    += rterm;
949:       ++k;
950:       /* Only need to evaluate every other point on refined levels */
951:       if (l != 1) ++k;
952:     } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */

954:     d1 = PetscLog10Real(PetscAbsReal(sum - osum));
955:     d2 = PetscLog10Real(PetscAbsReal(sum - psum));
956:     d3 = PetscLog10Real(maxTerm) - p;
957:     if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
958:     else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
959:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
960:   } while (d < digits && l < 12);
961:   *sol = sum;

963:   return(0);
964: }

966: #if defined(PETSC_HAVE_MPFR)
967: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
968: {
969:   const PetscInt  safetyFactor = 2;  /* Calculate abcissa until 2*p digits */
970:   PetscInt        l            = 0;  /* Level of refinement, h = 2^{-l} */
971:   mpfr_t          alpha;             /* Half-width of the integration interval */
972:   mpfr_t          beta;              /* Center of the integration interval */
973:   mpfr_t          h;                 /* Step size, length between x_k */
974:   mpfr_t          osum;              /* Integral on last level */
975:   mpfr_t          psum;              /* Integral on the level before the last level */
976:   mpfr_t          sum;               /* Integral on current level */
977:   mpfr_t          yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
978:   mpfr_t          lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
979:   mpfr_t          wk;                /* Quadrature weight at x_k */
980:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
981:   PetscInt        d;                 /* Digits of precision in the integral */
982:   mpfr_t          pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;

985:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
986:   /* Create high precision storage */
987:   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);
988:   /* Initialization */
989:   mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
990:   mpfr_set_d(beta,  0.5*(b+a), MPFR_RNDN);
991:   mpfr_set_d(osum,  0.0,       MPFR_RNDN);
992:   mpfr_set_d(psum,  0.0,       MPFR_RNDN);
993:   mpfr_set_d(h,     1.0,       MPFR_RNDN);
994:   mpfr_const_pi(pi2, MPFR_RNDN);
995:   mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
996:   /* Center term */
997:   func(0.5*(b+a), &lval);
998:   mpfr_set(sum, pi2, MPFR_RNDN);
999:   mpfr_mul(sum, sum, alpha, MPFR_RNDN);
1000:   mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
1001:   /* */
1002:   do {
1003:     PetscReal d1, d2, d3, d4;
1004:     PetscInt  k = 1;

1006:     ++l;
1007:     mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
1008:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1009:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1010:     mpfr_set(psum, osum, MPFR_RNDN);
1011:     mpfr_set(osum,  sum, MPFR_RNDN);
1012:     mpfr_mul_d(h,   h,   0.5, MPFR_RNDN);
1013:     mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
1014:     do {
1015:       mpfr_set_si(kh, k, MPFR_RNDN);
1016:       mpfr_mul(kh, kh, h, MPFR_RNDN);
1017:       /* Weight */
1018:       mpfr_set(wk, h, MPFR_RNDN);
1019:       mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
1020:       mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
1021:       mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
1022:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
1023:       mpfr_sqr(tmp, tmp, MPFR_RNDN);
1024:       mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
1025:       mpfr_div(wk, wk, tmp, MPFR_RNDN);
1026:       /* Abscissa */
1027:       mpfr_set_d(yk, 1.0, MPFR_RNDZ);
1028:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
1029:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1030:       mpfr_exp(tmp, msinh, MPFR_RNDN);
1031:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1032:       /* Quadrature points */
1033:       mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
1034:       mpfr_mul(lx, lx, alpha, MPFR_RNDU);
1035:       mpfr_add(lx, lx, beta, MPFR_RNDU);
1036:       mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
1037:       mpfr_mul(rx, rx, alpha, MPFR_RNDD);
1038:       mpfr_add(rx, rx, beta, MPFR_RNDD);
1039:       /* Evaluation */
1040:       func(mpfr_get_d(lx, MPFR_RNDU), &lval);
1041:       func(mpfr_get_d(rx, MPFR_RNDD), &rval);
1042:       /* Update */
1043:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1044:       mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
1045:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
1046:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1047:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1048:       mpfr_set(curTerm, tmp, MPFR_RNDN);
1049:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1050:       mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
1051:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
1052:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1053:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1054:       mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
1055:       ++k;
1056:       /* Only need to evaluate every other point on refined levels */
1057:       if (l != 1) ++k;
1058:       mpfr_log10(tmp, wk, MPFR_RNDN);
1059:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1060:     } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
1061:     mpfr_sub(tmp, sum, osum, MPFR_RNDN);
1062:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1063:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1064:     d1 = mpfr_get_d(tmp, MPFR_RNDN);
1065:     mpfr_sub(tmp, sum, psum, MPFR_RNDN);
1066:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1067:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1068:     d2 = mpfr_get_d(tmp, MPFR_RNDN);
1069:     mpfr_log10(tmp, maxTerm, MPFR_RNDN);
1070:     d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
1071:     mpfr_log10(tmp, curTerm, MPFR_RNDN);
1072:     d4 = mpfr_get_d(tmp, MPFR_RNDN);
1073:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1074:   } while (d < digits && l < 8);
1075:   *sol = mpfr_get_d(sum, MPFR_RNDN);
1076:   /* Cleanup */
1077:   mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
1078:   return(0);
1079: }
1080: #else

1082: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1083: {
1084:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
1085: }
1086: #endif

1088: /* Overwrites A. Can only handle full-rank problems with m>=n
1089:  * A in column-major format
1090:  * Ainv in row-major format
1091:  * tau has length m
1092:  * worksize must be >= max(1,n)
1093:  */
1094: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1095: {
1097:   PetscBLASInt   M,N,K,lda,ldb,ldwork,info;
1098:   PetscScalar    *A,*Ainv,*R,*Q,Alpha;

1101: #if defined(PETSC_USE_COMPLEX)
1102:   {
1103:     PetscInt i,j;
1104:     PetscMalloc2(m*n,&A,m*n,&Ainv);
1105:     for (j=0; j<n; j++) {
1106:       for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1107:     }
1108:     mstride = m;
1109:   }
1110: #else
1111:   A = A_in;
1112:   Ainv = Ainv_out;
1113: #endif

1115:   PetscBLASIntCast(m,&M);
1116:   PetscBLASIntCast(n,&N);
1117:   PetscBLASIntCast(mstride,&lda);
1118:   PetscBLASIntCast(worksize,&ldwork);
1119:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1120:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1121:   PetscFPTrapPop();
1122:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1123:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1125:   /* Extract an explicit representation of Q */
1126:   Q = Ainv;
1127:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1128:   K = N;                        /* full rank */
1129:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1130:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

1132:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1133:   Alpha = 1.0;
1134:   ldb = lda;
1135:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1136:   /* Ainv is Q, overwritten with inverse */

1138: #if defined(PETSC_USE_COMPLEX)
1139:   {
1140:     PetscInt i;
1141:     for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1142:     PetscFree2(A,Ainv);
1143:   }
1144: #endif
1145:   return(0);
1146: }

1148: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1149: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1150: {
1152:   PetscReal      *Bv;
1153:   PetscInt       i,j;

1156:   PetscMalloc1((ninterval+1)*ndegree,&Bv);
1157:   /* Point evaluation of L_p on all the source vertices */
1158:   PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1159:   /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1160:   for (i=0; i<ninterval; i++) {
1161:     for (j=0; j<ndegree; j++) {
1162:       if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1163:       else           B[i*ndegree+j]   = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1164:     }
1165:   }
1166:   PetscFree(Bv);
1167:   return(0);
1168: }

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

1173:    Not Collective

1175:    Input Arguments:
1176: +  degree - degree of reconstruction polynomial
1177: .  nsource - number of source intervals
1178: .  sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1179: .  ntarget - number of target intervals
1180: -  targetx - sorted coordinates of target cell boundaries (length ntarget+1)

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

1185:    Level: advanced

1187: .seealso: PetscDTLegendreEval()
1188: @*/
1189: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1190: {
1192:   PetscInt       i,j,k,*bdegrees,worksize;
1193:   PetscReal      xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1194:   PetscScalar    *tau,*work;

1200:   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);
1201: #if defined(PETSC_USE_DEBUG)
1202:   for (i=0; i<nsource; i++) {
1203:     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]);
1204:   }
1205:   for (i=0; i<ntarget; i++) {
1206:     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]);
1207:   }
1208: #endif
1209:   xmin = PetscMin(sourcex[0],targetx[0]);
1210:   xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1211:   center = (xmin + xmax)/2;
1212:   hscale = (xmax - xmin)/2;
1213:   worksize = nsource;
1214:   PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1215:   PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1216:   for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1217:   for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1218:   PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1219:   PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1220:   for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1221:   PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1222:   for (i=0; i<ntarget; i++) {
1223:     PetscReal rowsum = 0;
1224:     for (j=0; j<nsource; j++) {
1225:       PetscReal sum = 0;
1226:       for (k=0; k<degree+1; k++) {
1227:         sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1228:       }
1229:       R[i*nsource+j] = sum;
1230:       rowsum += sum;
1231:     }
1232:     for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1233:   }
1234:   PetscFree4(bdegrees,sourcey,Bsource,work);
1235:   PetscFree4(tau,Bsinv,targety,Btarget);
1236:   return(0);
1237: }