Actual source code: dt.c

petsc-3.9.2 2018-05-20
Report Typos and Errors
  1: /* Discretization tools */

  3: #include <petscconf.h>
  4: #if defined(PETSC_HAVE_MATHIMF_H)
  5: #include <mathimf.h>           /* this needs to be included before math.h */
  6: #endif
  7: #ifdef PETSC_HAVE_MPFR
  8: #include <mpfr.h>
  9: #endif

 11:  #include <petscdt.h>
 12:  #include <petscblaslapack.h>
 13:  #include <petsc/private/petscimpl.h>
 14:  #include <petsc/private/dtimpl.h>
 15:  #include <petscviewer.h>
 16:  #include <petscdmplex.h>
 17:  #include <petscdmshell.h>

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

 29: /*@
 30:   PetscQuadratureCreate - Create a PetscQuadrature object

 32:   Collective on MPI_Comm

 34:   Input Parameter:
 35: . comm - The communicator for the PetscQuadrature object

 37:   Output Parameter:
 38: . q  - The PetscQuadrature object

 40:   Level: beginner

 42: .keywords: PetscQuadrature, quadrature, create
 43: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
 44: @*/
 45: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
 46: {

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

 62: /*@
 63:   PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object

 65:   Collective on PetscQuadrature

 67:   Input Parameter:
 68: . q  - The PetscQuadrature object

 70:   Output Parameter:
 71: . r  - The new PetscQuadrature object

 73:   Level: beginner

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

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

 99: /*@
100:   PetscQuadratureDestroy - Destroys a PetscQuadrature object

102:   Collective on PetscQuadrature

104:   Input Parameter:
105: . q  - The PetscQuadrature object

107:   Level: beginner

109: .keywords: PetscQuadrature, quadrature, destroy
110: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
111: @*/
112: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
113: {

117:   if (!*q) return(0);
119:   if (--((PetscObject)(*q))->refct > 0) {
120:     *q = NULL;
121:     return(0);
122:   }
123:   PetscFree((*q)->points);
124:   PetscFree((*q)->weights);
125:   PetscHeaderDestroy(q);
126:   return(0);
127: }

129: /*@
130:   PetscQuadratureGetOrder - Return the order of the method

132:   Not collective

134:   Input Parameter:
135: . q - The PetscQuadrature object

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

140:   Level: intermediate

142: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
143: @*/
144: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
145: {
149:   *order = q->order;
150:   return(0);
151: }

153: /*@
154:   PetscQuadratureSetOrder - Return the order of the method

156:   Not collective

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

162:   Level: intermediate

164: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
165: @*/
166: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
167: {
170:   q->order = order;
171:   return(0);
172: }

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

177:   Not collective

179:   Input Parameter:
180: . q - The PetscQuadrature object

182:   Output Parameter:
183: . Nc - The number of components

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

187:   Level: intermediate

189: .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
190: @*/
191: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
192: {
196:   *Nc = q->Nc;
197:   return(0);
198: }

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

203:   Not collective

205:   Input Parameters:
206: + q  - The PetscQuadrature object
207: - Nc - The number of components

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

211:   Level: intermediate

213: .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
214: @*/
215: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
216: {
219:   q->Nc = Nc;
220:   return(0);
221: }

223: /*@C
224:   PetscQuadratureGetData - Returns the data defining the quadrature

226:   Not collective

228:   Input Parameter:
229: . q  - The PetscQuadrature object

231:   Output Parameters:
232: + dim - The spatial dimension
233: , Nc - The number of components
234: . npoints - The number of quadrature points
235: . points - The coordinates of each quadrature point
236: - weights - The weight of each quadrature point

238:   Level: intermediate

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

242: .keywords: PetscQuadrature, quadrature
243: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
244: @*/
245: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
246: {
249:   if (dim) {
251:     *dim = q->dim;
252:   }
253:   if (Nc) {
255:     *Nc = q->Nc;
256:   }
257:   if (npoints) {
259:     *npoints = q->numPoints;
260:   }
261:   if (points) {
263:     *points = q->points;
264:   }
265:   if (weights) {
267:     *weights = q->weights;
268:   }
269:   return(0);
270: }

272: /*@C
273:   PetscQuadratureSetData - Sets the data defining the quadrature

275:   Not collective

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

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

287:   Level: intermediate

289: .keywords: PetscQuadrature, quadrature
290: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
291: @*/
292: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
293: {
296:   if (dim >= 0)     q->dim       = dim;
297:   if (Nc >= 0)      q->Nc        = Nc;
298:   if (npoints >= 0) q->numPoints = npoints;
299:   if (points) {
301:     q->points = points;
302:   }
303:   if (weights) {
305:     q->weights = weights;
306:   }
307:   return(0);
308: }

310: /*@C
311:   PetscQuadratureView - Views a PetscQuadrature object

313:   Collective on PetscQuadrature

315:   Input Parameters:
316: + q  - The PetscQuadrature object
317: - viewer - The PetscViewer object

319:   Level: beginner

321: .keywords: PetscQuadrature, quadrature, view
322: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
323: @*/
324: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
325: {
326:   PetscInt       q, d, c;

330:   PetscObjectPrintClassNamePrefixType((PetscObject)quad,viewer);
331:   if (quad->Nc > 1) {PetscViewerASCIIPrintf(viewer, "Quadrature on %D points with %D components\n  (", quad->numPoints, quad->Nc);}
332:   else              {PetscViewerASCIIPrintf(viewer, "Quadrature on %D points\n  (", quad->numPoints);}
333:   for (q = 0; q < quad->numPoints; ++q) {
334:     for (d = 0; d < quad->dim; ++d) {
335:       if (d) PetscViewerASCIIPrintf(viewer, ", ");
336:       PetscViewerASCIIPrintf(viewer, "%g\n", (double)quad->points[q*quad->dim+d]);
337:     }
338:     if (quad->Nc > 1) {
339:       PetscViewerASCIIPrintf(viewer, ") (");
340:       for (c = 0; c < quad->Nc; ++c) {
341:         if (c) PetscViewerASCIIPrintf(viewer, ", ");
342:         PetscViewerASCIIPrintf(viewer, "%g", (double)quad->weights[q*quad->Nc+c]);
343:       }
344:       PetscViewerASCIIPrintf(viewer, ")\n");
345:     } else {
346:       PetscViewerASCIIPrintf(viewer, ") %g\n", (double)quad->weights[q]);
347:     }
348:   }
349:   return(0);
350: }

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

355:   Not collective

357:   Input Parameter:
358: + q - The original PetscQuadrature
359: . numSubelements - The number of subelements the original element is divided into
360: . v0 - An array of the initial points for each subelement
361: - jac - An array of the Jacobian mappings from the reference to each subelement

363:   Output Parameters:
364: . dim - The dimension

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

368:  Not available from Fortran

370:   Level: intermediate

372: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
373: @*/
374: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
375: {
376:   const PetscReal *points,    *weights;
377:   PetscReal       *pointsRef, *weightsRef;
378:   PetscInt         dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
379:   PetscErrorCode   ierr;

386:   PetscQuadratureCreate(PETSC_COMM_SELF, qref);
387:   PetscQuadratureGetOrder(q, &order);
388:   PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
389:   npointsRef = npoints*numSubelements;
390:   PetscMalloc1(npointsRef*dim,&pointsRef);
391:   PetscMalloc1(npointsRef*Nc, &weightsRef);
392:   for (c = 0; c < numSubelements; ++c) {
393:     for (p = 0; p < npoints; ++p) {
394:       for (d = 0; d < dim; ++d) {
395:         pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
396:         for (e = 0; e < dim; ++e) {
397:           pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
398:         }
399:       }
400:       /* Could also use detJ here */
401:       for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
402:     }
403:   }
404:   PetscQuadratureSetOrder(*qref, order);
405:   PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
406:   return(0);
407: }

409: /*@
410:    PetscDTLegendreEval - evaluate Legendre polynomial at points

412:    Not Collective

414:    Input Arguments:
415: +  npoints - number of spatial points to evaluate at
416: .  points - array of locations to evaluate at
417: .  ndegree - number of basis degrees to evaluate
418: -  degrees - sorted array of degrees to evaluate

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

425:    Level: intermediate

427: .seealso: PetscDTGaussQuadrature()
428: @*/
429: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
430: {
431:   PetscInt i,maxdegree;

434:   if (!npoints || !ndegree) return(0);
435:   maxdegree = degrees[ndegree-1];
436:   for (i=0; i<npoints; i++) {
437:     PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
438:     PetscInt  j,k;
439:     x    = points[i];
440:     pm2  = 0;
441:     pm1  = 1;
442:     pd2  = 0;
443:     pd1  = 0;
444:     pdd2 = 0;
445:     pdd1 = 0;
446:     k    = 0;
447:     if (degrees[k] == 0) {
448:       if (B) B[i*ndegree+k] = pm1;
449:       if (D) D[i*ndegree+k] = pd1;
450:       if (D2) D2[i*ndegree+k] = pdd1;
451:       k++;
452:     }
453:     for (j=1; j<=maxdegree; j++,k++) {
454:       PetscReal p,d,dd;
455:       p    = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
456:       d    = pd2 + (2*j-1)*pm1;
457:       dd   = pdd2 + (2*j-1)*pd1;
458:       pm2  = pm1;
459:       pm1  = p;
460:       pd2  = pd1;
461:       pd1  = d;
462:       pdd2 = pdd1;
463:       pdd1 = dd;
464:       if (degrees[k] == j) {
465:         if (B) B[i*ndegree+k] = p;
466:         if (D) D[i*ndegree+k] = d;
467:         if (D2) D2[i*ndegree+k] = dd;
468:       }
469:     }
470:   }
471:   return(0);
472: }

474: /*@
475:    PetscDTGaussQuadrature - create Gauss quadrature

477:    Not Collective

479:    Input Arguments:
480: +  npoints - number of points
481: .  a - left end of interval (often-1)
482: -  b - right end of interval (often +1)

484:    Output Arguments:
485: +  x - quadrature points
486: -  w - quadrature weights

488:    Level: intermediate

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

493: .seealso: PetscDTLegendreEval()
494: @*/
495: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
496: {
498:   PetscInt       i;
499:   PetscReal      *work;
500:   PetscScalar    *Z;
501:   PetscBLASInt   N,LDZ,info;

504:   PetscCitationsRegister(GaussCitation, &GaussCite);
505:   /* Set up the Golub-Welsch system */
506:   for (i=0; i<npoints; i++) {
507:     x[i] = 0;                   /* diagonal is 0 */
508:     if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
509:   }
510:   PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
511:   PetscBLASIntCast(npoints,&N);
512:   LDZ  = N;
513:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
514:   PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
515:   PetscFPTrapPop();
516:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");

518:   for (i=0; i<(npoints+1)/2; i++) {
519:     PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
520:     x[i]           = (a+b)/2 - y*(b-a)/2;
521:     if (x[i] == -0.0) x[i] = 0.0;
522:     x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;

524:     w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
525:   }
526:   PetscFree2(Z,work);
527:   return(0);
528: }

530: /*@
531:   PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature

533:   Not Collective

535:   Input Arguments:
536: + dim     - The spatial dimension
537: . Nc      - The number of components
538: . npoints - number of points in one dimension
539: . a       - left end of interval (often-1)
540: - b       - right end of interval (often +1)

542:   Output Argument:
543: . q - A PetscQuadrature object

545:   Level: intermediate

547: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
548: @*/
549: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
550: {
551:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
552:   PetscReal     *x, *w, *xw, *ww;

556:   PetscMalloc1(totpoints*dim,&x);
557:   PetscMalloc1(totpoints*Nc,&w);
558:   /* Set up the Golub-Welsch system */
559:   switch (dim) {
560:   case 0:
561:     PetscFree(x);
562:     PetscFree(w);
563:     PetscMalloc1(1, &x);
564:     PetscMalloc1(Nc, &w);
565:     x[0] = 0.0;
566:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
567:     break;
568:   case 1:
569:     PetscMalloc1(npoints,&ww);
570:     PetscDTGaussQuadrature(npoints, a, b, x, ww);
571:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
572:     PetscFree(ww);
573:     break;
574:   case 2:
575:     PetscMalloc2(npoints,&xw,npoints,&ww);
576:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
577:     for (i = 0; i < npoints; ++i) {
578:       for (j = 0; j < npoints; ++j) {
579:         x[(i*npoints+j)*dim+0] = xw[i];
580:         x[(i*npoints+j)*dim+1] = xw[j];
581:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
582:       }
583:     }
584:     PetscFree2(xw,ww);
585:     break;
586:   case 3:
587:     PetscMalloc2(npoints,&xw,npoints,&ww);
588:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
589:     for (i = 0; i < npoints; ++i) {
590:       for (j = 0; j < npoints; ++j) {
591:         for (k = 0; k < npoints; ++k) {
592:           x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
593:           x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
594:           x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
595:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
596:         }
597:       }
598:     }
599:     PetscFree2(xw,ww);
600:     break;
601:   default:
602:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
603:   }
604:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
605:   PetscQuadratureSetOrder(*q, npoints-1);
606:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
607:   return(0);
608: }

610: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
611:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
612: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
613: {
614:   PetscReal f = 1.0;
615:   PetscInt  i;

618:   for (i = 1; i < n+1; ++i) f *= i;
619:   *factorial = f;
620:   return(0);
621: }

623: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
624:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
625: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
626: {
627:   PetscReal apb, pn1, pn2;
628:   PetscInt  k;

631:   if (!n) {*P = 1.0; return(0);}
632:   if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
633:   apb = a + b;
634:   pn2 = 1.0;
635:   pn1 = 0.5 * (a - b + (apb + 2.0) * x);
636:   *P  = 0.0;
637:   for (k = 2; k < n+1; ++k) {
638:     PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
639:     PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
640:     PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
641:     PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);

643:     a2  = a2 / a1;
644:     a3  = a3 / a1;
645:     a4  = a4 / a1;
646:     *P  = (a2 + a3 * x) * pn1 - a4 * pn2;
647:     pn2 = pn1;
648:     pn1 = *P;
649:   }
650:   return(0);
651: }

653: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
654: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
655: {
656:   PetscReal      nP;

660:   if (!n) {*P = 0.0; return(0);}
661:   PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
662:   *P   = 0.5 * (a + b + n + 1) * nP;
663:   return(0);
664: }

666: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
667: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
668: {
670:   *xi  = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
671:   *eta = y;
672:   return(0);
673: }

675: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
676: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
677: {
679:   *xi   = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
680:   *eta  = 0.5  * (1.0 + y) * (1.0 - z) - 1.0;
681:   *zeta = z;
682:   return(0);
683: }

685: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
686: {
687:   PetscInt       maxIter = 100;
688:   PetscReal      eps     = 1.0e-8;
689:   PetscReal      a1, a2, a3, a4, a5, a6;
690:   PetscInt       k;


695:   a1      = PetscPowReal(2.0, a+b+1);
696: #if defined(PETSC_HAVE_TGAMMA)
697:   a2      = PetscTGamma(a + npoints + 1);
698:   a3      = PetscTGamma(b + npoints + 1);
699:   a4      = PetscTGamma(a + b + npoints + 1);
700: #else
701:   {
702:     PetscInt ia, ib;

704:     ia = (PetscInt) a;
705:     ib = (PetscInt) b;
706:     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 */
707:       PetscDTFactorial_Internal(ia + npoints, &a2);
708:       PetscDTFactorial_Internal(ib + npoints, &a3);
709:       PetscDTFactorial_Internal(ia + ib + npoints, &a4);
710:     } else {
711:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
712:     }
713:   }
714: #endif

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

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

729:       for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
730:       PetscDTComputeJacobi(a, b, npoints, r, &f);
731:       PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
732:       delta = f / (fp - f * s);
733:       r     = r - delta;
734:       if (PetscAbsReal(delta) < eps) break;
735:     }
736:     x[k] = r;
737:     PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
738:     w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
739:   }
740:   return(0);
741: }

743: /*@
744:   PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex

746:   Not Collective

748:   Input Arguments:
749: + dim     - The simplex dimension
750: . Nc      - The number of components
751: . npoints - The number of points in one dimension
752: . a       - left end of interval (often-1)
753: - b       - right end of interval (often +1)

755:   Output Argument:
756: . q - A PetscQuadrature object

758:   Level: intermediate

760:   References:
761: .  1. - Karniadakis and Sherwin.  FIAT

763: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
764: @*/
765: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
766: {
767:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints;
768:   PetscReal     *px, *wx, *py, *wy, *pz, *wz, *x, *w;
769:   PetscInt       i, j, k, c;

773:   if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
774:   PetscMalloc1(totpoints*dim, &x);
775:   PetscMalloc1(totpoints*Nc, &w);
776:   switch (dim) {
777:   case 0:
778:     PetscFree(x);
779:     PetscFree(w);
780:     PetscMalloc1(1, &x);
781:     PetscMalloc1(Nc, &w);
782:     x[0] = 0.0;
783:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
784:     break;
785:   case 1:
786:     PetscMalloc1(npoints,&wx);
787:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, wx);
788:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i];
789:     PetscFree(wx);
790:     break;
791:   case 2:
792:     PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);
793:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
794:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
795:     for (i = 0; i < npoints; ++i) {
796:       for (j = 0; j < npoints; ++j) {
797:         PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);
798:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j];
799:       }
800:     }
801:     PetscFree4(px,wx,py,wy);
802:     break;
803:   case 3:
804:     PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);
805:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
806:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
807:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);
808:     for (i = 0; i < npoints; ++i) {
809:       for (j = 0; j < npoints; ++j) {
810:         for (k = 0; k < npoints; ++k) {
811:           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]);
812:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k];
813:         }
814:       }
815:     }
816:     PetscFree6(px,wx,py,wy,pz,wz);
817:     break;
818:   default:
819:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
820:   }
821:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
822:   PetscQuadratureSetOrder(*q, npoints-1);
823:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
824:   return(0);
825: }

827: /*@
828:   PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell

830:   Not Collective

832:   Input Arguments:
833: + dim   - The cell dimension
834: . level - The number of points in one dimension, 2^l
835: . a     - left end of interval (often-1)
836: - b     - right end of interval (often +1)

838:   Output Argument:
839: . q - A PetscQuadrature object

841:   Level: intermediate

843: .seealso: PetscDTGaussTensorQuadrature()
844: @*/
845: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
846: {
847:   const PetscInt  p     = 16;                        /* Digits of precision in the evaluation */
848:   const PetscReal alpha = (b-a)/2.;                  /* Half-width of the integration interval */
849:   const PetscReal beta  = (b+a)/2.;                  /* Center of the integration interval */
850:   const PetscReal h     = PetscPowReal(2.0, -level); /* Step size, length between x_k */
851:   PetscReal       xk;                                /* Quadrature point x_k on reference domain [-1, 1] */
852:   PetscReal       wk    = 0.5*PETSC_PI;              /* Quadrature weight at x_k */
853:   PetscReal      *x, *w;
854:   PetscInt        K, k, npoints;
855:   PetscErrorCode  ierr;

858:   if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
859:   if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
860:   /* Find K such that the weights are < 32 digits of precision */
861:   for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
862:     wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
863:   }
864:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
865:   PetscQuadratureSetOrder(*q, 2*K+1);
866:   npoints = 2*K-1;
867:   PetscMalloc1(npoints*dim, &x);
868:   PetscMalloc1(npoints, &w);
869:   /* Center term */
870:   x[0] = beta;
871:   w[0] = 0.5*alpha*PETSC_PI;
872:   for (k = 1; k < K; ++k) {
873:     wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
874:     xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
875:     x[2*k-1] = -alpha*xk+beta;
876:     w[2*k-1] = wk;
877:     x[2*k+0] =  alpha*xk+beta;
878:     w[2*k+0] = wk;
879:   }
880:   PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
881:   return(0);
882: }

884: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
885: {
886:   const PetscInt  p     = 16;        /* Digits of precision in the evaluation */
887:   const PetscReal alpha = (b-a)/2.;  /* Half-width of the integration interval */
888:   const PetscReal beta  = (b+a)/2.;  /* Center of the integration interval */
889:   PetscReal       h     = 1.0;       /* Step size, length between x_k */
890:   PetscInt        l     = 0;         /* Level of refinement, h = 2^{-l} */
891:   PetscReal       osum  = 0.0;       /* Integral on last level */
892:   PetscReal       psum  = 0.0;       /* Integral on the level before the last level */
893:   PetscReal       sum;               /* Integral on current level */
894:   PetscReal       yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
895:   PetscReal       lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
896:   PetscReal       wk;                /* Quadrature weight at x_k */
897:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
898:   PetscInt        d;                 /* Digits of precision in the integral */

901:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
902:   /* Center term */
903:   func(beta, &lval);
904:   sum = 0.5*alpha*PETSC_PI*lval;
905:   /* */
906:   do {
907:     PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
908:     PetscInt  k = 1;

910:     ++l;
911:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
912:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
913:     psum = osum;
914:     osum = sum;
915:     h   *= 0.5;
916:     sum *= 0.5;
917:     do {
918:       wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
919:       yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
920:       lx = -alpha*(1.0 - yk)+beta;
921:       rx =  alpha*(1.0 - yk)+beta;
922:       func(lx, &lval);
923:       func(rx, &rval);
924:       lterm   = alpha*wk*lval;
925:       maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
926:       sum    += lterm;
927:       rterm   = alpha*wk*rval;
928:       maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
929:       sum    += rterm;
930:       ++k;
931:       /* Only need to evaluate every other point on refined levels */
932:       if (l != 1) ++k;
933:     } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */

935:     d1 = PetscLog10Real(PetscAbsReal(sum - osum));
936:     d2 = PetscLog10Real(PetscAbsReal(sum - psum));
937:     d3 = PetscLog10Real(maxTerm) - p;
938:     if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
939:     else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
940:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
941:   } while (d < digits && l < 12);
942:   *sol = sum;

944:   return(0);
945: }

947: #ifdef PETSC_HAVE_MPFR
948: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
949: {
950:   const PetscInt  safetyFactor = 2;  /* Calculate abcissa until 2*p digits */
951:   PetscInt        l            = 0;  /* Level of refinement, h = 2^{-l} */
952:   mpfr_t          alpha;             /* Half-width of the integration interval */
953:   mpfr_t          beta;              /* Center of the integration interval */
954:   mpfr_t          h;                 /* Step size, length between x_k */
955:   mpfr_t          osum;              /* Integral on last level */
956:   mpfr_t          psum;              /* Integral on the level before the last level */
957:   mpfr_t          sum;               /* Integral on current level */
958:   mpfr_t          yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
959:   mpfr_t          lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
960:   mpfr_t          wk;                /* Quadrature weight at x_k */
961:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
962:   PetscInt        d;                 /* Digits of precision in the integral */
963:   mpfr_t          pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;

966:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
967:   /* Create high precision storage */
968:   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);
969:   /* Initialization */
970:   mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
971:   mpfr_set_d(beta,  0.5*(b+a), MPFR_RNDN);
972:   mpfr_set_d(osum,  0.0,       MPFR_RNDN);
973:   mpfr_set_d(psum,  0.0,       MPFR_RNDN);
974:   mpfr_set_d(h,     1.0,       MPFR_RNDN);
975:   mpfr_const_pi(pi2, MPFR_RNDN);
976:   mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
977:   /* Center term */
978:   func(0.5*(b+a), &lval);
979:   mpfr_set(sum, pi2, MPFR_RNDN);
980:   mpfr_mul(sum, sum, alpha, MPFR_RNDN);
981:   mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
982:   /* */
983:   do {
984:     PetscReal d1, d2, d3, d4;
985:     PetscInt  k = 1;

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

1063: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1064: {
1065:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
1066: }
1067: #endif

1069: /* Overwrites A. Can only handle full-rank problems with m>=n
1070:  * A in column-major format
1071:  * Ainv in row-major format
1072:  * tau has length m
1073:  * worksize must be >= max(1,n)
1074:  */
1075: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1076: {
1078:   PetscBLASInt   M,N,K,lda,ldb,ldwork,info;
1079:   PetscScalar    *A,*Ainv,*R,*Q,Alpha;

1082: #if defined(PETSC_USE_COMPLEX)
1083:   {
1084:     PetscInt i,j;
1085:     PetscMalloc2(m*n,&A,m*n,&Ainv);
1086:     for (j=0; j<n; j++) {
1087:       for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1088:     }
1089:     mstride = m;
1090:   }
1091: #else
1092:   A = A_in;
1093:   Ainv = Ainv_out;
1094: #endif

1096:   PetscBLASIntCast(m,&M);
1097:   PetscBLASIntCast(n,&N);
1098:   PetscBLASIntCast(mstride,&lda);
1099:   PetscBLASIntCast(worksize,&ldwork);
1100:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1101:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1102:   PetscFPTrapPop();
1103:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1104:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1106:   /* Extract an explicit representation of Q */
1107:   Q = Ainv;
1108:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1109:   K = N;                        /* full rank */
1110:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1111:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

1113:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1114:   Alpha = 1.0;
1115:   ldb = lda;
1116:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1117:   /* Ainv is Q, overwritten with inverse */

1119: #if defined(PETSC_USE_COMPLEX)
1120:   {
1121:     PetscInt i;
1122:     for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1123:     PetscFree2(A,Ainv);
1124:   }
1125: #endif
1126:   return(0);
1127: }

1129: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1130: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1131: {
1133:   PetscReal      *Bv;
1134:   PetscInt       i,j;

1137:   PetscMalloc1((ninterval+1)*ndegree,&Bv);
1138:   /* Point evaluation of L_p on all the source vertices */
1139:   PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1140:   /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1141:   for (i=0; i<ninterval; i++) {
1142:     for (j=0; j<ndegree; j++) {
1143:       if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1144:       else           B[i*ndegree+j]   = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1145:     }
1146:   }
1147:   PetscFree(Bv);
1148:   return(0);
1149: }

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

1154:    Not Collective

1156:    Input Arguments:
1157: +  degree - degree of reconstruction polynomial
1158: .  nsource - number of source intervals
1159: .  sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1160: .  ntarget - number of target intervals
1161: -  targetx - sorted coordinates of target cell boundaries (length ntarget+1)

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

1166:    Level: advanced

1168: .seealso: PetscDTLegendreEval()
1169: @*/
1170: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1171: {
1173:   PetscInt       i,j,k,*bdegrees,worksize;
1174:   PetscReal      xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1175:   PetscScalar    *tau,*work;

1181:   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);
1182: #if defined(PETSC_USE_DEBUG)
1183:   for (i=0; i<nsource; i++) {
1184:     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]);
1185:   }
1186:   for (i=0; i<ntarget; i++) {
1187:     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]);
1188:   }
1189: #endif
1190:   xmin = PetscMin(sourcex[0],targetx[0]);
1191:   xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1192:   center = (xmin + xmax)/2;
1193:   hscale = (xmax - xmin)/2;
1194:   worksize = nsource;
1195:   PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1196:   PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1197:   for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1198:   for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1199:   PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1200:   PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1201:   for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1202:   PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1203:   for (i=0; i<ntarget; i++) {
1204:     PetscReal rowsum = 0;
1205:     for (j=0; j<nsource; j++) {
1206:       PetscReal sum = 0;
1207:       for (k=0; k<degree+1; k++) {
1208:         sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1209:       }
1210:       R[i*nsource+j] = sum;
1211:       rowsum += sum;
1212:     }
1213:     for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1214:   }
1215:   PetscFree4(bdegrees,sourcey,Bsource,work);
1216:   PetscFree4(tau,Bsinv,targety,Btarget);
1217:   return(0);
1218: }