Actual source code: dt.c

petsc-master 2017-04-29
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)->order     = -1;
 55:   (*q)->numPoints = 0;
 56:   (*q)->points    = NULL;
 57:   (*q)->weights   = NULL;
 58:   return(0);
 59: }

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

 64:   Collective on PetscQuadrature

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

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

 72:   Level: beginner

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

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

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

101:   Collective on PetscQuadrature

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

106:   Level: beginner

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

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

128: /*@
129:   PetscQuadratureGetOrder - Return the quadrature information

131:   Not collective

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

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

139:   Output Parameter:

141:   Level: intermediate

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

154: /*@
155:   PetscQuadratureSetOrder - Return the quadrature information

157:   Not collective

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

163:   Level: intermediate

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

175: /*@C
176:   PetscQuadratureGetData - Returns the data defining the quadrature

178:   Not collective

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

183:   Output Parameters:
184: + dim - The spatial dimension
185: . npoints - The number of quadrature points
186: . points - The coordinates of each quadrature point
187: - weights - The weight of each quadrature point

189:   Level: intermediate

191: .keywords: PetscQuadrature, quadrature
192: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
193: @*/
194: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
195: {
198:   if (dim) {
200:     *dim = q->dim;
201:   }
202:   if (npoints) {
204:     *npoints = q->numPoints;
205:   }
206:   if (points) {
208:     *points = q->points;
209:   }
210:   if (weights) {
212:     *weights = q->weights;
213:   }
214:   return(0);
215: }

217: /*@C
218:   PetscQuadratureSetData - Sets the data defining the quadrature

220:   Not collective

222:   Input Parameters:
223: + q  - The PetscQuadrature object
224: . dim - The spatial dimension
225: . npoints - The number of quadrature points
226: . points - The coordinates of each quadrature point
227: - weights - The weight of each quadrature point

229:   Level: intermediate

231: .keywords: PetscQuadrature, quadrature
232: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
233: @*/
234: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
235: {
238:   if (dim >= 0)     q->dim       = dim;
239:   if (npoints >= 0) q->numPoints = npoints;
240:   if (points) {
242:     q->points = points;
243:   }
244:   if (weights) {
246:     q->weights = weights;
247:   }
248:   return(0);
249: }

251: /*@C
252:   PetscQuadratureView - Views a PetscQuadrature object

254:   Collective on PetscQuadrature

256:   Input Parameters:
257: + q  - The PetscQuadrature object
258: - viewer - The PetscViewer object

260:   Level: beginner

262: .keywords: PetscQuadrature, quadrature, view
263: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
264: @*/
265: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
266: {
267:   PetscInt       q, d;

271:   PetscObjectPrintClassNamePrefixType((PetscObject)quad,viewer);
272:   PetscViewerASCIIPrintf(viewer, "Quadrature on %d points\n  (", quad->numPoints);
273:   for (q = 0; q < quad->numPoints; ++q) {
274:     for (d = 0; d < quad->dim; ++d) {
275:       if (d) PetscViewerASCIIPrintf(viewer, ", ");
276:       PetscViewerASCIIPrintf(viewer, "%g\n", (double)quad->points[q*quad->dim+d]);
277:     }
278:     PetscViewerASCIIPrintf(viewer, ") %g\n", (double)quad->weights[q]);
279:   }
280:   return(0);
281: }

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

286:   Not collective

288:   Input Parameter:
289: + q - The original PetscQuadrature
290: . numSubelements - The number of subelements the original element is divided into
291: . v0 - An array of the initial points for each subelement
292: - jac - An array of the Jacobian mappings from the reference to each subelement

294:   Output Parameters:
295: . dim - The dimension

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

299:   Level: intermediate

301: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
302: @*/
303: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
304: {
305:   const PetscReal *points,    *weights;
306:   PetscReal       *pointsRef, *weightsRef;
307:   PetscInt         dim, order, npoints, npointsRef, c, p, d, e;
308:   PetscErrorCode   ierr;

315:   PetscQuadratureCreate(PETSC_COMM_SELF, qref);
316:   PetscQuadratureGetOrder(q, &order);
317:   PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
318:   npointsRef = npoints*numSubelements;
319:   PetscMalloc1(npointsRef*dim,&pointsRef);
320:   PetscMalloc1(npointsRef,&weightsRef);
321:   for (c = 0; c < numSubelements; ++c) {
322:     for (p = 0; p < npoints; ++p) {
323:       for (d = 0; d < dim; ++d) {
324:         pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
325:         for (e = 0; e < dim; ++e) {
326:           pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
327:         }
328:       }
329:       /* Could also use detJ here */
330:       weightsRef[c*npoints+p] = weights[p]/numSubelements;
331:     }
332:   }
333:   PetscQuadratureSetOrder(*qref, order);
334:   PetscQuadratureSetData(*qref, dim, npointsRef, pointsRef, weightsRef);
335:   return(0);
336: }

338: /*@
339:    PetscDTLegendreEval - evaluate Legendre polynomial at points

341:    Not Collective

343:    Input Arguments:
344: +  npoints - number of spatial points to evaluate at
345: .  points - array of locations to evaluate at
346: .  ndegree - number of basis degrees to evaluate
347: -  degrees - sorted array of degrees to evaluate

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

354:    Level: intermediate

356: .seealso: PetscDTGaussQuadrature()
357: @*/
358: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
359: {
360:   PetscInt i,maxdegree;

363:   if (!npoints || !ndegree) return(0);
364:   maxdegree = degrees[ndegree-1];
365:   for (i=0; i<npoints; i++) {
366:     PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
367:     PetscInt  j,k;
368:     x    = points[i];
369:     pm2  = 0;
370:     pm1  = 1;
371:     pd2  = 0;
372:     pd1  = 0;
373:     pdd2 = 0;
374:     pdd1 = 0;
375:     k    = 0;
376:     if (degrees[k] == 0) {
377:       if (B) B[i*ndegree+k] = pm1;
378:       if (D) D[i*ndegree+k] = pd1;
379:       if (D2) D2[i*ndegree+k] = pdd1;
380:       k++;
381:     }
382:     for (j=1; j<=maxdegree; j++,k++) {
383:       PetscReal p,d,dd;
384:       p    = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
385:       d    = pd2 + (2*j-1)*pm1;
386:       dd   = pdd2 + (2*j-1)*pd1;
387:       pm2  = pm1;
388:       pm1  = p;
389:       pd2  = pd1;
390:       pd1  = d;
391:       pdd2 = pdd1;
392:       pdd1 = dd;
393:       if (degrees[k] == j) {
394:         if (B) B[i*ndegree+k] = p;
395:         if (D) D[i*ndegree+k] = d;
396:         if (D2) D2[i*ndegree+k] = dd;
397:       }
398:     }
399:   }
400:   return(0);
401: }

403: /*@
404:    PetscDTGaussQuadrature - create Gauss quadrature

406:    Not Collective

408:    Input Arguments:
409: +  npoints - number of points
410: .  a - left end of interval (often-1)
411: -  b - right end of interval (often +1)

413:    Output Arguments:
414: +  x - quadrature points
415: -  w - quadrature weights

417:    Level: intermediate

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

422: .seealso: PetscDTLegendreEval()
423: @*/
424: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
425: {
427:   PetscInt       i;
428:   PetscReal      *work;
429:   PetscScalar    *Z;
430:   PetscBLASInt   N,LDZ,info;

433:   PetscCitationsRegister(GaussCitation, &GaussCite);
434:   /* Set up the Golub-Welsch system */
435:   for (i=0; i<npoints; i++) {
436:     x[i] = 0;                   /* diagonal is 0 */
437:     if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
438:   }
439:   PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
440:   PetscBLASIntCast(npoints,&N);
441:   LDZ  = N;
442:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
443:   PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
444:   PetscFPTrapPop();
445:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");

447:   for (i=0; i<(npoints+1)/2; i++) {
448:     PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
449:     x[i]           = (a+b)/2 - y*(b-a)/2;
450:     if (x[i] == -0.0) x[i] = 0.0;
451:     x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;

453:     w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
454:   }
455:   PetscFree2(Z,work);
456:   return(0);
457: }

459: /*@
460:   PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature

462:   Not Collective

464:   Input Arguments:
465: + dim     - The spatial dimension
466: . npoints - number of points in one dimension
467: . a       - left end of interval (often-1)
468: - b       - right end of interval (often +1)

470:   Output Argument:
471: . q - A PetscQuadrature object

473:   Level: intermediate

475: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
476: @*/
477: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
478: {
479:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k;
480:   PetscReal     *x, *w, *xw, *ww;

484:   PetscMalloc1(totpoints*dim,&x);
485:   PetscMalloc1(totpoints,&w);
486:   /* Set up the Golub-Welsch system */
487:   switch (dim) {
488:   case 0:
489:     PetscFree(x);
490:     PetscFree(w);
491:     PetscMalloc1(1, &x);
492:     PetscMalloc1(1, &w);
493:     x[0] = 0.0;
494:     w[0] = 1.0;
495:     break;
496:   case 1:
497:     PetscDTGaussQuadrature(npoints, a, b, x, w);
498:     break;
499:   case 2:
500:     PetscMalloc2(npoints,&xw,npoints,&ww);
501:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
502:     for (i = 0; i < npoints; ++i) {
503:       for (j = 0; j < npoints; ++j) {
504:         x[(i*npoints+j)*dim+0] = xw[i];
505:         x[(i*npoints+j)*dim+1] = xw[j];
506:         w[i*npoints+j]         = ww[i] * ww[j];
507:       }
508:     }
509:     PetscFree2(xw,ww);
510:     break;
511:   case 3:
512:     PetscMalloc2(npoints,&xw,npoints,&ww);
513:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
514:     for (i = 0; i < npoints; ++i) {
515:       for (j = 0; j < npoints; ++j) {
516:         for (k = 0; k < npoints; ++k) {
517:           x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
518:           x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
519:           x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
520:           w[(i*npoints+j)*npoints+k]         = ww[i] * ww[j] * ww[k];
521:         }
522:       }
523:     }
524:     PetscFree2(xw,ww);
525:     break;
526:   default:
527:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
528:   }
529:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
530:   PetscQuadratureSetOrder(*q, npoints);
531:   PetscQuadratureSetData(*q, dim, totpoints, x, w);
532:   return(0);
533: }

535: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
536:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
537: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
538: {
539:   PetscReal f = 1.0;
540:   PetscInt  i;

543:   for (i = 1; i < n+1; ++i) f *= i;
544:   *factorial = f;
545:   return(0);
546: }

548: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
549:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
550: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
551: {
552:   PetscReal apb, pn1, pn2;
553:   PetscInt  k;

556:   if (!n) {*P = 1.0; return(0);}
557:   if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
558:   apb = a + b;
559:   pn2 = 1.0;
560:   pn1 = 0.5 * (a - b + (apb + 2.0) * x);
561:   *P  = 0.0;
562:   for (k = 2; k < n+1; ++k) {
563:     PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
564:     PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
565:     PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
566:     PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);

568:     a2  = a2 / a1;
569:     a3  = a3 / a1;
570:     a4  = a4 / a1;
571:     *P  = (a2 + a3 * x) * pn1 - a4 * pn2;
572:     pn2 = pn1;
573:     pn1 = *P;
574:   }
575:   return(0);
576: }

578: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
579: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
580: {
581:   PetscReal      nP;

585:   if (!n) {*P = 0.0; return(0);}
586:   PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
587:   *P   = 0.5 * (a + b + n + 1) * nP;
588:   return(0);
589: }

591: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
592: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
593: {
595:   *xi  = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
596:   *eta = y;
597:   return(0);
598: }

600: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
601: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
602: {
604:   *xi   = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
605:   *eta  = 0.5  * (1.0 + y) * (1.0 - z) - 1.0;
606:   *zeta = z;
607:   return(0);
608: }

610: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
611: {
612:   PetscInt       maxIter = 100;
613:   PetscReal      eps     = 1.0e-8;
614:   PetscReal      a1, a2, a3, a4, a5, a6;
615:   PetscInt       k;


620:   a1      = PetscPowReal(2.0, a+b+1);
621: #if defined(PETSC_HAVE_TGAMMA)
622:   a2      = PetscTGamma(a + npoints + 1);
623:   a3      = PetscTGamma(b + npoints + 1);
624:   a4      = PetscTGamma(a + b + npoints + 1);
625: #else
626:   {
627:     PetscInt ia, ib;

629:     ia = (PetscInt) a;
630:     ib = (PetscInt) b;
631:     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 */
632:       PetscDTFactorial_Internal(ia + npoints, &a2);
633:       PetscDTFactorial_Internal(ib + npoints, &a3);
634:       PetscDTFactorial_Internal(ia + ib + npoints, &a4);
635:     } else {
636:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
637:     }
638:   }
639: #endif

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

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

654:       for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
655:       PetscDTComputeJacobi(a, b, npoints, r, &f);
656:       PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
657:       delta = f / (fp - f * s);
658:       r     = r - delta;
659:       if (PetscAbsReal(delta) < eps) break;
660:     }
661:     x[k] = r;
662:     PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
663:     w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
664:   }
665:   return(0);
666: }

668: /*@C
669:   PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex

671:   Not Collective

673:   Input Arguments:
674: + dim   - The simplex dimension
675: . order - The number of points in one dimension
676: . a     - left end of interval (often-1)
677: - b     - right end of interval (often +1)

679:   Output Argument:
680: . q - A PetscQuadrature object

682:   Level: intermediate

684:   References:
685: .  1. - Karniadakis and Sherwin.  FIAT

687: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
688: @*/
689: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt order, PetscReal a, PetscReal b, PetscQuadrature *q)
690: {
691:   PetscInt       npoints = dim > 1 ? dim > 2 ? order*PetscSqr(order) : PetscSqr(order) : order;
692:   PetscReal     *px, *wx, *py, *wy, *pz, *wz, *x, *w;
693:   PetscInt       i, j, k;

697:   if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
698:   PetscMalloc1(npoints*dim, &x);
699:   PetscMalloc1(npoints, &w);
700:   switch (dim) {
701:   case 0:
702:     PetscFree(x);
703:     PetscFree(w);
704:     PetscMalloc1(1, &x);
705:     PetscMalloc1(1, &w);
706:     x[0] = 0.0;
707:     w[0] = 1.0;
708:     break;
709:   case 1:
710:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, x, w);
711:     break;
712:   case 2:
713:     PetscMalloc4(order,&px,order,&wx,order,&py,order,&wy);
714:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
715:     PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
716:     for (i = 0; i < order; ++i) {
717:       for (j = 0; j < order; ++j) {
718:         PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*order+j)*2+0], &x[(i*order+j)*2+1]);
719:         w[i*order+j] = 0.5 * wx[i] * wy[j];
720:       }
721:     }
722:     PetscFree4(px,wx,py,wy);
723:     break;
724:   case 3:
725:     PetscMalloc6(order,&px,order,&wx,order,&py,order,&wy,order,&pz,order,&wz);
726:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
727:     PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
728:     PetscDTGaussJacobiQuadrature1D_Internal(order, 2.0, 0.0, pz, wz);
729:     for (i = 0; i < order; ++i) {
730:       for (j = 0; j < order; ++j) {
731:         for (k = 0; k < order; ++k) {
732:           PetscDTMapCubeToTetrahedron_Internal(px[i], py[j], pz[k], &x[((i*order+j)*order+k)*3+0], &x[((i*order+j)*order+k)*3+1], &x[((i*order+j)*order+k)*3+2]);
733:           w[(i*order+j)*order+k] = 0.125 * wx[i] * wy[j] * wz[k];
734:         }
735:       }
736:     }
737:     PetscFree6(px,wx,py,wy,pz,wz);
738:     break;
739:   default:
740:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
741:   }
742:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
743:   PetscQuadratureSetOrder(*q, order);
744:   PetscQuadratureSetData(*q, dim, npoints, x, w);
745:   return(0);
746: }

748: /*@C
749:   PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell

751:   Not Collective

753:   Input Arguments:
754: + dim   - The cell dimension
755: . level - The number of points in one dimension, 2^l
756: . a     - left end of interval (often-1)
757: - b     - right end of interval (often +1)

759:   Output Argument:
760: . q - A PetscQuadrature object

762:   Level: intermediate

764: .seealso: PetscDTGaussTensorQuadrature()
765: @*/
766: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
767: {
768:   const PetscInt  p     = 16;                        /* Digits of precision in the evaluation */
769:   const PetscReal alpha = (b-a)/2.;                  /* Half-width of the integration interval */
770:   const PetscReal beta  = (b+a)/2.;                  /* Center of the integration interval */
771:   const PetscReal h     = PetscPowReal(2.0, -level); /* Step size, length between x_k */
772:   PetscReal       xk;                                /* Quadrature point x_k on reference domain [-1, 1] */
773:   PetscReal       wk    = 0.5*PETSC_PI;              /* Quadrature weight at x_k */
774:   PetscReal      *x, *w;
775:   PetscInt        K, k, npoints;
776:   PetscErrorCode  ierr;

779:   if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
780:   if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
781:   /* Find K such that the weights are < 32 digits of precision */
782:   for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
783:     wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
784:   }
785:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
786:   PetscQuadratureSetOrder(*q, 2*K+1);
787:   npoints = 2*K-1;
788:   PetscMalloc1(npoints*dim, &x);
789:   PetscMalloc1(npoints, &w);
790:   /* Center term */
791:   x[0] = beta;
792:   w[0] = 0.5*alpha*PETSC_PI;
793:   for (k = 1; k < K; ++k) {
794:     wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
795:     xk = tanh(0.5*PETSC_PI*PetscSinhReal(k*h));
796:     x[2*k-1] = -alpha*xk+beta;
797:     w[2*k-1] = wk;
798:     x[2*k+0] =  alpha*xk+beta;
799:     w[2*k+0] = wk;
800:   }
801:   PetscQuadratureSetData(*q, dim, npoints, x, w);
802:   return(0);
803: }

805: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
806: {
807:   const PetscInt  p     = 16;        /* Digits of precision in the evaluation */
808:   const PetscReal alpha = (b-a)/2.;  /* Half-width of the integration interval */
809:   const PetscReal beta  = (b+a)/2.;  /* Center of the integration interval */
810:   PetscReal       h     = 1.0;       /* Step size, length between x_k */
811:   PetscInt        l     = 0;         /* Level of refinement, h = 2^{-l} */
812:   PetscReal       osum  = 0.0;       /* Integral on last level */
813:   PetscReal       psum  = 0.0;       /* Integral on the level before the last level */
814:   PetscReal       sum;               /* Integral on current level */
815:   PetscReal       yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
816:   PetscReal       lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
817:   PetscReal       wk;                /* Quadrature weight at x_k */
818:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
819:   PetscInt        d;                 /* Digits of precision in the integral */

822:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
823:   /* Center term */
824:   func(beta, &lval);
825:   sum = 0.5*alpha*PETSC_PI*lval;
826:   /* */
827:   do {
828:     PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
829:     PetscInt  k = 1;

831:     ++l;
832:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
833:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
834:     psum = osum;
835:     osum = sum;
836:     h   *= 0.5;
837:     sum *= 0.5;
838:     do {
839:       wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
840:       yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
841:       lx = -alpha*(1.0 - yk)+beta;
842:       rx =  alpha*(1.0 - yk)+beta;
843:       func(lx, &lval);
844:       func(rx, &rval);
845:       lterm   = alpha*wk*lval;
846:       maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
847:       sum    += lterm;
848:       rterm   = alpha*wk*rval;
849:       maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
850:       sum    += rterm;
851:       ++k;
852:       /* Only need to evaluate every other point on refined levels */
853:       if (l != 1) ++k;
854:     } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */

856:     d1 = PetscLog10Real(PetscAbsReal(sum - osum));
857:     d2 = PetscLog10Real(PetscAbsReal(sum - psum));
858:     d3 = PetscLog10Real(maxTerm) - p;
859:     if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
860:     else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
861:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
862:   } while (d < digits && l < 12);
863:   *sol = sum;

865:   return(0);
866: }

868: #ifdef PETSC_HAVE_MPFR
869: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
870: {
871:   const PetscInt  safetyFactor = 2;  /* Calculate abcissa until 2*p digits */
872:   PetscInt        l            = 0;  /* Level of refinement, h = 2^{-l} */
873:   mpfr_t          alpha;             /* Half-width of the integration interval */
874:   mpfr_t          beta;              /* Center of the integration interval */
875:   mpfr_t          h;                 /* Step size, length between x_k */
876:   mpfr_t          osum;              /* Integral on last level */
877:   mpfr_t          psum;              /* Integral on the level before the last level */
878:   mpfr_t          sum;               /* Integral on current level */
879:   mpfr_t          yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
880:   mpfr_t          lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
881:   mpfr_t          wk;                /* Quadrature weight at x_k */
882:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
883:   PetscInt        d;                 /* Digits of precision in the integral */
884:   mpfr_t          pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;

887:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
888:   /* Create high precision storage */
889:   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);
890:   /* Initialization */
891:   mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
892:   mpfr_set_d(beta,  0.5*(b+a), MPFR_RNDN);
893:   mpfr_set_d(osum,  0.0,       MPFR_RNDN);
894:   mpfr_set_d(psum,  0.0,       MPFR_RNDN);
895:   mpfr_set_d(h,     1.0,       MPFR_RNDN);
896:   mpfr_const_pi(pi2, MPFR_RNDN);
897:   mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
898:   /* Center term */
899:   func(0.5*(b+a), &lval);
900:   mpfr_set(sum, pi2, MPFR_RNDN);
901:   mpfr_mul(sum, sum, alpha, MPFR_RNDN);
902:   mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
903:   /* */
904:   do {
905:     PetscReal d1, d2, d3, d4;
906:     PetscInt  k = 1;

908:     ++l;
909:     mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
910:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
911:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
912:     mpfr_set(psum, osum, MPFR_RNDN);
913:     mpfr_set(osum,  sum, MPFR_RNDN);
914:     mpfr_mul_d(h,   h,   0.5, MPFR_RNDN);
915:     mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
916:     do {
917:       mpfr_set_si(kh, k, MPFR_RNDN);
918:       mpfr_mul(kh, kh, h, MPFR_RNDN);
919:       /* Weight */
920:       mpfr_set(wk, h, MPFR_RNDN);
921:       mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
922:       mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
923:       mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
924:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
925:       mpfr_sqr(tmp, tmp, MPFR_RNDN);
926:       mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
927:       mpfr_div(wk, wk, tmp, MPFR_RNDN);
928:       /* Abscissa */
929:       mpfr_set_d(yk, 1.0, MPFR_RNDZ);
930:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
931:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
932:       mpfr_exp(tmp, msinh, MPFR_RNDN);
933:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
934:       /* Quadrature points */
935:       mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
936:       mpfr_mul(lx, lx, alpha, MPFR_RNDU);
937:       mpfr_add(lx, lx, beta, MPFR_RNDU);
938:       mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
939:       mpfr_mul(rx, rx, alpha, MPFR_RNDD);
940:       mpfr_add(rx, rx, beta, MPFR_RNDD);
941:       /* Evaluation */
942:       func(mpfr_get_d(lx, MPFR_RNDU), &lval);
943:       func(mpfr_get_d(rx, MPFR_RNDD), &rval);
944:       /* Update */
945:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
946:       mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
947:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
948:       mpfr_abs(tmp, tmp, MPFR_RNDN);
949:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
950:       mpfr_set(curTerm, tmp, MPFR_RNDN);
951:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
952:       mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
953:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
954:       mpfr_abs(tmp, tmp, MPFR_RNDN);
955:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
956:       mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
957:       ++k;
958:       /* Only need to evaluate every other point on refined levels */
959:       if (l != 1) ++k;
960:       mpfr_log10(tmp, wk, MPFR_RNDN);
961:       mpfr_abs(tmp, tmp, MPFR_RNDN);
962:     } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
963:     mpfr_sub(tmp, sum, osum, MPFR_RNDN);
964:     mpfr_abs(tmp, tmp, MPFR_RNDN);
965:     mpfr_log10(tmp, tmp, MPFR_RNDN);
966:     d1 = mpfr_get_d(tmp, MPFR_RNDN);
967:     mpfr_sub(tmp, sum, psum, MPFR_RNDN);
968:     mpfr_abs(tmp, tmp, MPFR_RNDN);
969:     mpfr_log10(tmp, tmp, MPFR_RNDN);
970:     d2 = mpfr_get_d(tmp, MPFR_RNDN);
971:     mpfr_log10(tmp, maxTerm, MPFR_RNDN);
972:     d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
973:     mpfr_log10(tmp, curTerm, MPFR_RNDN);
974:     d4 = mpfr_get_d(tmp, MPFR_RNDN);
975:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
976:   } while (d < digits && l < 8);
977:   *sol = mpfr_get_d(sum, MPFR_RNDN);
978:   /* Cleanup */
979:   mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
980:   return(0);
981: }
982: #else

984: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
985: {
986:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
987: }
988: #endif

990: /* Overwrites A. Can only handle full-rank problems with m>=n
991:  * A in column-major format
992:  * Ainv in row-major format
993:  * tau has length m
994:  * worksize must be >= max(1,n)
995:  */
996: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
997: {
999:   PetscBLASInt   M,N,K,lda,ldb,ldwork,info;
1000:   PetscScalar    *A,*Ainv,*R,*Q,Alpha;

1003: #if defined(PETSC_USE_COMPLEX)
1004:   {
1005:     PetscInt i,j;
1006:     PetscMalloc2(m*n,&A,m*n,&Ainv);
1007:     for (j=0; j<n; j++) {
1008:       for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1009:     }
1010:     mstride = m;
1011:   }
1012: #else
1013:   A = A_in;
1014:   Ainv = Ainv_out;
1015: #endif

1017:   PetscBLASIntCast(m,&M);
1018:   PetscBLASIntCast(n,&N);
1019:   PetscBLASIntCast(mstride,&lda);
1020:   PetscBLASIntCast(worksize,&ldwork);
1021:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1022:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1023:   PetscFPTrapPop();
1024:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1025:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1027:   /* Extract an explicit representation of Q */
1028:   Q = Ainv;
1029:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1030:   K = N;                        /* full rank */
1031:   PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1032:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

1034:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1035:   Alpha = 1.0;
1036:   ldb = lda;
1037:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1038:   /* Ainv is Q, overwritten with inverse */

1040: #if defined(PETSC_USE_COMPLEX)
1041:   {
1042:     PetscInt i;
1043:     for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1044:     PetscFree2(A,Ainv);
1045:   }
1046: #endif
1047:   return(0);
1048: }

1050: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1051: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1052: {
1054:   PetscReal      *Bv;
1055:   PetscInt       i,j;

1058:   PetscMalloc1((ninterval+1)*ndegree,&Bv);
1059:   /* Point evaluation of L_p on all the source vertices */
1060:   PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1061:   /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1062:   for (i=0; i<ninterval; i++) {
1063:     for (j=0; j<ndegree; j++) {
1064:       if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1065:       else           B[i*ndegree+j]   = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1066:     }
1067:   }
1068:   PetscFree(Bv);
1069:   return(0);
1070: }

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

1075:    Not Collective

1077:    Input Arguments:
1078: +  degree - degree of reconstruction polynomial
1079: .  nsource - number of source intervals
1080: .  sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1081: .  ntarget - number of target intervals
1082: -  targetx - sorted coordinates of target cell boundaries (length ntarget+1)

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

1087:    Level: advanced

1089: .seealso: PetscDTLegendreEval()
1090: @*/
1091: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1092: {
1094:   PetscInt       i,j,k,*bdegrees,worksize;
1095:   PetscReal      xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1096:   PetscScalar    *tau,*work;

1102:   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);
1103: #if defined(PETSC_USE_DEBUG)
1104:   for (i=0; i<nsource; i++) {
1105:     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]);
1106:   }
1107:   for (i=0; i<ntarget; i++) {
1108:     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]);
1109:   }
1110: #endif
1111:   xmin = PetscMin(sourcex[0],targetx[0]);
1112:   xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1113:   center = (xmin + xmax)/2;
1114:   hscale = (xmax - xmin)/2;
1115:   worksize = nsource;
1116:   PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1117:   PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1118:   for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1119:   for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1120:   PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1121:   PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1122:   for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1123:   PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1124:   for (i=0; i<ntarget; i++) {
1125:     PetscReal rowsum = 0;
1126:     for (j=0; j<nsource; j++) {
1127:       PetscReal sum = 0;
1128:       for (k=0; k<degree+1; k++) {
1129:         sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1130:       }
1131:       R[i*nsource+j] = sum;
1132:       rowsum += sum;
1133:     }
1134:     for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1135:   }
1136:   PetscFree4(bdegrees,sourcey,Bsource,work);
1137:   PetscFree4(tau,Bsinv,targety,Btarget);
1138:   return(0);
1139: }