Actual source code: dt.c

petsc-master 2016-12-06
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";

 31: /*@
 32:   PetscQuadratureCreate - Create a PetscQuadrature object

 34:   Collective on MPI_Comm

 36:   Input Parameter:
 37: . comm - The communicator for the PetscQuadrature object

 39:   Output Parameter:
 40: . q  - The PetscQuadrature object

 42:   Level: beginner

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

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

 65: /*@
 66:   PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object

 68:   Collective on PetscQuadrature

 70:   Input Parameter:
 71: . q  - The PetscQuadrature object

 73:   Output Parameter:
 74: . r  - The new PetscQuadrature object

 76:   Level: beginner

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

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

104: /*@
105:   PetscQuadratureDestroy - Destroys a PetscQuadrature object

107:   Collective on PetscQuadrature

109:   Input Parameter:
110: . q  - The PetscQuadrature object

112:   Level: beginner

114: .keywords: PetscQuadrature, quadrature, destroy
115: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
116: @*/
117: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
118: {

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

136: /*@
137:   PetscQuadratureGetOrder - Return the quadrature information

139:   Not collective

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

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

147:   Output Parameter:

149:   Level: intermediate

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

164: /*@
165:   PetscQuadratureSetOrder - Return the quadrature information

167:   Not collective

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

173:   Level: intermediate

175: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
176: @*/
177: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
178: {
181:   q->order = order;
182:   return(0);
183: }

187: /*@C
188:   PetscQuadratureGetData - Returns the data defining the quadrature

190:   Not collective

192:   Input Parameter:
193: . q  - The PetscQuadrature object

195:   Output Parameters:
196: + dim - The spatial dimension
197: . npoints - The number of quadrature points
198: . points - The coordinates of each quadrature point
199: - weights - The weight of each quadrature point

201:   Level: intermediate

203: .keywords: PetscQuadrature, quadrature
204: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
205: @*/
206: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
207: {
210:   if (dim) {
212:     *dim = q->dim;
213:   }
214:   if (npoints) {
216:     *npoints = q->numPoints;
217:   }
218:   if (points) {
220:     *points = q->points;
221:   }
222:   if (weights) {
224:     *weights = q->weights;
225:   }
226:   return(0);
227: }

231: /*@C
232:   PetscQuadratureSetData - Sets the data defining the quadrature

234:   Not collective

236:   Input Parameters:
237: + q  - The PetscQuadrature object
238: . dim - The spatial dimension
239: . npoints - The number of quadrature points
240: . points - The coordinates of each quadrature point
241: - weights - The weight of each quadrature point

243:   Level: intermediate

245: .keywords: PetscQuadrature, quadrature
246: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
247: @*/
248: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
249: {
252:   if (dim >= 0)     q->dim       = dim;
253:   if (npoints >= 0) q->numPoints = npoints;
254:   if (points) {
256:     q->points = points;
257:   }
258:   if (weights) {
260:     q->weights = weights;
261:   }
262:   return(0);
263: }

267: /*@C
268:   PetscQuadratureView - Views a PetscQuadrature object

270:   Collective on PetscQuadrature

272:   Input Parameters:
273: + q  - The PetscQuadrature object
274: - viewer - The PetscViewer object

276:   Level: beginner

278: .keywords: PetscQuadrature, quadrature, view
279: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
280: @*/
281: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
282: {
283:   PetscInt       q, d;

287:   PetscObjectPrintClassNamePrefixType((PetscObject)quad,viewer);
288:   PetscViewerASCIIPrintf(viewer, "Quadrature on %d points\n  (", quad->numPoints);
289:   for (q = 0; q < quad->numPoints; ++q) {
290:     for (d = 0; d < quad->dim; ++d) {
291:       if (d) PetscViewerASCIIPrintf(viewer, ", ");
292:       PetscViewerASCIIPrintf(viewer, "%g\n", (double)quad->points[q*quad->dim+d]);
293:     }
294:     PetscViewerASCIIPrintf(viewer, ") %g\n", (double)quad->weights[q]);
295:   }
296:   return(0);
297: }

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

304:   Not collective

306:   Input Parameter:
307: + q - The original PetscQuadrature
308: . numSubelements - The number of subelements the original element is divided into
309: . v0 - An array of the initial points for each subelement
310: - jac - An array of the Jacobian mappings from the reference to each subelement

312:   Output Parameters:
313: . dim - The dimension

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

317:   Level: intermediate

319: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
320: @*/
321: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
322: {
323:   const PetscReal *points,    *weights;
324:   PetscReal       *pointsRef, *weightsRef;
325:   PetscInt         dim, order, npoints, npointsRef, c, p, d, e;
326:   PetscErrorCode   ierr;

333:   PetscQuadratureCreate(PETSC_COMM_SELF, qref);
334:   PetscQuadratureGetOrder(q, &order);
335:   PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
336:   npointsRef = npoints*numSubelements;
337:   PetscMalloc1(npointsRef*dim,&pointsRef);
338:   PetscMalloc1(npointsRef,&weightsRef);
339:   for (c = 0; c < numSubelements; ++c) {
340:     for (p = 0; p < npoints; ++p) {
341:       for (d = 0; d < dim; ++d) {
342:         pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
343:         for (e = 0; e < dim; ++e) {
344:           pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
345:         }
346:       }
347:       /* Could also use detJ here */
348:       weightsRef[c*npoints+p] = weights[p]/numSubelements;
349:     }
350:   }
351:   PetscQuadratureSetOrder(*qref, order);
352:   PetscQuadratureSetData(*qref, dim, npointsRef, pointsRef, weightsRef);
353:   return(0);
354: }

358: /*@
359:    PetscDTLegendreEval - evaluate Legendre polynomial at points

361:    Not Collective

363:    Input Arguments:
364: +  npoints - number of spatial points to evaluate at
365: .  points - array of locations to evaluate at
366: .  ndegree - number of basis degrees to evaluate
367: -  degrees - sorted array of degrees to evaluate

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

374:    Level: intermediate

376: .seealso: PetscDTGaussQuadrature()
377: @*/
378: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
379: {
380:   PetscInt i,maxdegree;

383:   if (!npoints || !ndegree) return(0);
384:   maxdegree = degrees[ndegree-1];
385:   for (i=0; i<npoints; i++) {
386:     PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
387:     PetscInt  j,k;
388:     x    = points[i];
389:     pm2  = 0;
390:     pm1  = 1;
391:     pd2  = 0;
392:     pd1  = 0;
393:     pdd2 = 0;
394:     pdd1 = 0;
395:     k    = 0;
396:     if (degrees[k] == 0) {
397:       if (B) B[i*ndegree+k] = pm1;
398:       if (D) D[i*ndegree+k] = pd1;
399:       if (D2) D2[i*ndegree+k] = pdd1;
400:       k++;
401:     }
402:     for (j=1; j<=maxdegree; j++,k++) {
403:       PetscReal p,d,dd;
404:       p    = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
405:       d    = pd2 + (2*j-1)*pm1;
406:       dd   = pdd2 + (2*j-1)*pd1;
407:       pm2  = pm1;
408:       pm1  = p;
409:       pd2  = pd1;
410:       pd1  = d;
411:       pdd2 = pdd1;
412:       pdd1 = dd;
413:       if (degrees[k] == j) {
414:         if (B) B[i*ndegree+k] = p;
415:         if (D) D[i*ndegree+k] = d;
416:         if (D2) D2[i*ndegree+k] = dd;
417:       }
418:     }
419:   }
420:   return(0);
421: }

425: /*@
426:    PetscDTGaussQuadrature - create Gauss quadrature

428:    Not Collective

430:    Input Arguments:
431: +  npoints - number of points
432: .  a - left end of interval (often-1)
433: -  b - right end of interval (often +1)

435:    Output Arguments:
436: +  x - quadrature points
437: -  w - quadrature weights

439:    Level: intermediate

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

444: .seealso: PetscDTLegendreEval()
445: @*/
446: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
447: {
449:   PetscInt       i;
450:   PetscReal      *work;
451:   PetscScalar    *Z;
452:   PetscBLASInt   N,LDZ,info;

455:   PetscCitationsRegister(GaussCitation, &GaussCite);
456:   /* Set up the Golub-Welsch system */
457:   for (i=0; i<npoints; i++) {
458:     x[i] = 0;                   /* diagonal is 0 */
459:     if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
460:   }
461:   PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
462:   PetscBLASIntCast(npoints,&N);
463:   LDZ  = N;
464:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
465:   PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
466:   PetscFPTrapPop();
467:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");

469:   for (i=0; i<(npoints+1)/2; i++) {
470:     PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
471:     x[i]           = (a+b)/2 - y*(b-a)/2;
472:     if (x[i] == -0.0) x[i] = 0.0;
473:     x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;

475:     w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
476:   }
477:   PetscFree2(Z,work);
478:   return(0);
479: }

483: /*@
484:   PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature

486:   Not Collective

488:   Input Arguments:
489: + dim     - The spatial dimension
490: . npoints - number of points in one dimension
491: . a       - left end of interval (often-1)
492: - b       - right end of interval (often +1)

494:   Output Argument:
495: . q - A PetscQuadrature object

497:   Level: intermediate

499: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
500: @*/
501: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
502: {
503:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k;
504:   PetscReal     *x, *w, *xw, *ww;

508:   PetscMalloc1(totpoints*dim,&x);
509:   PetscMalloc1(totpoints,&w);
510:   /* Set up the Golub-Welsch system */
511:   switch (dim) {
512:   case 0:
513:     PetscFree(x);
514:     PetscFree(w);
515:     PetscMalloc1(1, &x);
516:     PetscMalloc1(1, &w);
517:     x[0] = 0.0;
518:     w[0] = 1.0;
519:     break;
520:   case 1:
521:     PetscDTGaussQuadrature(npoints, a, b, x, w);
522:     break;
523:   case 2:
524:     PetscMalloc2(npoints,&xw,npoints,&ww);
525:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
526:     for (i = 0; i < npoints; ++i) {
527:       for (j = 0; j < npoints; ++j) {
528:         x[(i*npoints+j)*dim+0] = xw[i];
529:         x[(i*npoints+j)*dim+1] = xw[j];
530:         w[i*npoints+j]         = ww[i] * ww[j];
531:       }
532:     }
533:     PetscFree2(xw,ww);
534:     break;
535:   case 3:
536:     PetscMalloc2(npoints,&xw,npoints,&ww);
537:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
538:     for (i = 0; i < npoints; ++i) {
539:       for (j = 0; j < npoints; ++j) {
540:         for (k = 0; k < npoints; ++k) {
541:           x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
542:           x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
543:           x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
544:           w[(i*npoints+j)*npoints+k]         = ww[i] * ww[j] * ww[k];
545:         }
546:       }
547:     }
548:     PetscFree2(xw,ww);
549:     break;
550:   default:
551:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
552:   }
553:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
554:   PetscQuadratureSetOrder(*q, npoints);
555:   PetscQuadratureSetData(*q, dim, totpoints, x, w);
556:   return(0);
557: }

561: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
562:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
563: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
564: {
565:   PetscReal f = 1.0;
566:   PetscInt  i;

569:   for (i = 1; i < n+1; ++i) f *= i;
570:   *factorial = f;
571:   return(0);
572: }

576: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
577:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
578: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
579: {
580:   PetscReal apb, pn1, pn2;
581:   PetscInt  k;

584:   if (!n) {*P = 1.0; return(0);}
585:   if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
586:   apb = a + b;
587:   pn2 = 1.0;
588:   pn1 = 0.5 * (a - b + (apb + 2.0) * x);
589:   *P  = 0.0;
590:   for (k = 2; k < n+1; ++k) {
591:     PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
592:     PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
593:     PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
594:     PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);

596:     a2  = a2 / a1;
597:     a3  = a3 / a1;
598:     a4  = a4 / a1;
599:     *P  = (a2 + a3 * x) * pn1 - a4 * pn2;
600:     pn2 = pn1;
601:     pn1 = *P;
602:   }
603:   return(0);
604: }

608: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
609: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
610: {
611:   PetscReal      nP;

615:   if (!n) {*P = 0.0; return(0);}
616:   PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
617:   *P   = 0.5 * (a + b + n + 1) * nP;
618:   return(0);
619: }

623: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
624: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
625: {
627:   *xi  = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
628:   *eta = y;
629:   return(0);
630: }

634: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
635: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
636: {
638:   *xi   = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
639:   *eta  = 0.5  * (1.0 + y) * (1.0 - z) - 1.0;
640:   *zeta = z;
641:   return(0);
642: }

646: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
647: {
648:   PetscInt       maxIter = 100;
649:   PetscReal      eps     = 1.0e-8;
650:   PetscReal      a1, a2, a3, a4, a5, a6;
651:   PetscInt       k;


656:   a1      = PetscPowReal(2.0, a+b+1);
657: #if defined(PETSC_HAVE_TGAMMA)
658:   a2      = PetscTGamma(a + npoints + 1);
659:   a3      = PetscTGamma(b + npoints + 1);
660:   a4      = PetscTGamma(a + b + npoints + 1);
661: #else
662:   {
663:     PetscInt ia, ib;

665:     ia = (PetscInt) a;
666:     ib = (PetscInt) b;
667:     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 */
668:       PetscDTFactorial_Internal(ia + npoints, &a2);
669:       PetscDTFactorial_Internal(ib + npoints, &a3);
670:       PetscDTFactorial_Internal(ia + ib + npoints, &a4);
671:     } else {
672:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
673:     }
674:   }
675: #endif

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

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

690:       for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
691:       PetscDTComputeJacobi(a, b, npoints, r, &f);
692:       PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
693:       delta = f / (fp - f * s);
694:       r     = r - delta;
695:       if (PetscAbsReal(delta) < eps) break;
696:     }
697:     x[k] = r;
698:     PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
699:     w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
700:   }
701:   return(0);
702: }

706: /*@C
707:   PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex

709:   Not Collective

711:   Input Arguments:
712: + dim   - The simplex dimension
713: . order - The number of points in one dimension
714: . a     - left end of interval (often-1)
715: - b     - right end of interval (often +1)

717:   Output Argument:
718: . q - A PetscQuadrature object

720:   Level: intermediate

722:   References:
723: .  1. - Karniadakis and Sherwin.  FIAT

725: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
726: @*/
727: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt order, PetscReal a, PetscReal b, PetscQuadrature *q)
728: {
729:   PetscInt       npoints = dim > 1 ? dim > 2 ? order*PetscSqr(order) : PetscSqr(order) : order;
730:   PetscReal     *px, *wx, *py, *wy, *pz, *wz, *x, *w;
731:   PetscInt       i, j, k;

735:   if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
736:   PetscMalloc1(npoints*dim, &x);
737:   PetscMalloc1(npoints, &w);
738:   switch (dim) {
739:   case 0:
740:     PetscFree(x);
741:     PetscFree(w);
742:     PetscMalloc1(1, &x);
743:     PetscMalloc1(1, &w);
744:     x[0] = 0.0;
745:     w[0] = 1.0;
746:     break;
747:   case 1:
748:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, x, w);
749:     break;
750:   case 2:
751:     PetscMalloc4(order,&px,order,&wx,order,&py,order,&wy);
752:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
753:     PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
754:     for (i = 0; i < order; ++i) {
755:       for (j = 0; j < order; ++j) {
756:         PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*order+j)*2+0], &x[(i*order+j)*2+1]);
757:         w[i*order+j] = 0.5 * wx[i] * wy[j];
758:       }
759:     }
760:     PetscFree4(px,wx,py,wy);
761:     break;
762:   case 3:
763:     PetscMalloc6(order,&px,order,&wx,order,&py,order,&wy,order,&pz,order,&wz);
764:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
765:     PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
766:     PetscDTGaussJacobiQuadrature1D_Internal(order, 2.0, 0.0, pz, wz);
767:     for (i = 0; i < order; ++i) {
768:       for (j = 0; j < order; ++j) {
769:         for (k = 0; k < order; ++k) {
770:           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]);
771:           w[(i*order+j)*order+k] = 0.125 * wx[i] * wy[j] * wz[k];
772:         }
773:       }
774:     }
775:     PetscFree6(px,wx,py,wy,pz,wz);
776:     break;
777:   default:
778:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
779:   }
780:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
781:   PetscQuadratureSetOrder(*q, order);
782:   PetscQuadratureSetData(*q, dim, npoints, x, w);
783:   return(0);
784: }

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

791:   Not Collective

793:   Input Arguments:
794: + dim   - The cell dimension
795: . level - The number of points in one dimension, 2^l
796: . a     - left end of interval (often-1)
797: - b     - right end of interval (often +1)

799:   Output Argument:
800: . q - A PetscQuadrature object

802:   Level: intermediate

804: .seealso: PetscDTGaussTensorQuadrature()
805: @*/
806: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
807: {
808:   const PetscInt  p     = 16;                        /* Digits of precision in the evaluation */
809:   const PetscReal alpha = (b-a)/2.;                  /* Half-width of the integration interval */
810:   const PetscReal beta  = (b+a)/2.;                  /* Center of the integration interval */
811:   const PetscReal h     = PetscPowReal(2.0, -level); /* Step size, length between x_k */
812:   PetscReal       xk;                                /* Quadrature point x_k on reference domain [-1, 1] */
813:   PetscReal       wk    = 0.5*PETSC_PI;              /* Quadrature weight at x_k */
814:   PetscReal      *x, *w;
815:   PetscInt        K, k, npoints;
816:   PetscErrorCode  ierr;

819:   if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
820:   if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
821:   /* Find K such that the weights are < 32 digits of precision */
822:   for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
823:     wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
824:   }
825:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
826:   PetscQuadratureSetOrder(*q, 2*K+1);
827:   npoints = 2*K-1;
828:   PetscMalloc1(npoints*dim, &x);
829:   PetscMalloc1(npoints, &w);
830:   /* Center term */
831:   x[0] = beta;
832:   w[0] = 0.5*alpha*PETSC_PI;
833:   for (k = 1; k < K; ++k) {
834:     wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
835:     xk = tanh(0.5*PETSC_PI*PetscSinhReal(k*h));
836:     x[2*k-1] = -alpha*xk+beta;
837:     w[2*k-1] = wk;
838:     x[2*k+0] =  alpha*xk+beta;
839:     w[2*k+0] = wk;
840:   }
841:   PetscQuadratureSetData(*q, dim, npoints, x, w);
842:   return(0);
843: }

847: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
848: {
849:   const PetscInt  p     = 16;        /* Digits of precision in the evaluation */
850:   const PetscReal alpha = (b-a)/2.;  /* Half-width of the integration interval */
851:   const PetscReal beta  = (b+a)/2.;  /* Center of the integration interval */
852:   PetscReal       h     = 1.0;       /* Step size, length between x_k */
853:   PetscInt        l     = 0;         /* Level of refinement, h = 2^{-l} */
854:   PetscReal       osum  = 0.0;       /* Integral on last level */
855:   PetscReal       psum  = 0.0;       /* Integral on the level before the last level */
856:   PetscReal       sum;               /* Integral on current level */
857:   PetscReal       yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
858:   PetscReal       lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
859:   PetscReal       wk;                /* Quadrature weight at x_k */
860:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
861:   PetscInt        d;                 /* Digits of precision in the integral */

864:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
865:   /* Center term */
866:   func(beta, &lval);
867:   sum = 0.5*alpha*PETSC_PI*lval;
868:   /* */
869:   do {
870:     PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
871:     PetscInt  k = 1;

873:     ++l;
874:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
875:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
876:     psum = osum;
877:     osum = sum;
878:     h   *= 0.5;
879:     sum *= 0.5;
880:     do {
881:       wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
882:       yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
883:       lx = -alpha*(1.0 - yk)+beta;
884:       rx =  alpha*(1.0 - yk)+beta;
885:       func(lx, &lval);
886:       func(rx, &rval);
887:       lterm   = alpha*wk*lval;
888:       maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
889:       sum    += lterm;
890:       rterm   = alpha*wk*rval;
891:       maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
892:       sum    += rterm;
893:       ++k;
894:       /* Only need to evaluate every other point on refined levels */
895:       if (l != 1) ++k;
896:     } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */

898:     d1 = PetscLog10Real(PetscAbsReal(sum - osum));
899:     d2 = PetscLog10Real(PetscAbsReal(sum - psum));
900:     d3 = PetscLog10Real(maxTerm) - p;
901:     if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
902:     else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
903:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
904:   } while (d < digits && l < 12);
905:   *sol = sum;

907:   return(0);
908: }

910: #ifdef PETSC_HAVE_MPFR
913: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
914: {
915:   const PetscInt  safetyFactor = 2;  /* Calculate abcissa until 2*p digits */
916:   PetscInt        l            = 0;  /* Level of refinement, h = 2^{-l} */
917:   mpfr_t          alpha;             /* Half-width of the integration interval */
918:   mpfr_t          beta;              /* Center of the integration interval */
919:   mpfr_t          h;                 /* Step size, length between x_k */
920:   mpfr_t          osum;              /* Integral on last level */
921:   mpfr_t          psum;              /* Integral on the level before the last level */
922:   mpfr_t          sum;               /* Integral on current level */
923:   mpfr_t          yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
924:   mpfr_t          lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
925:   mpfr_t          wk;                /* Quadrature weight at x_k */
926:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
927:   PetscInt        d;                 /* Digits of precision in the integral */
928:   mpfr_t          pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;

931:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
932:   /* Create high precision storage */
933:   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);
934:   /* Initialization */
935:   mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
936:   mpfr_set_d(beta,  0.5*(b+a), MPFR_RNDN);
937:   mpfr_set_d(osum,  0.0,       MPFR_RNDN);
938:   mpfr_set_d(psum,  0.0,       MPFR_RNDN);
939:   mpfr_set_d(h,     1.0,       MPFR_RNDN);
940:   mpfr_const_pi(pi2, MPFR_RNDN);
941:   mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
942:   /* Center term */
943:   func(0.5*(b+a), &lval);
944:   mpfr_set(sum, pi2, MPFR_RNDN);
945:   mpfr_mul(sum, sum, alpha, MPFR_RNDN);
946:   mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
947:   /* */
948:   do {
949:     PetscReal d1, d2, d3, d4;
950:     PetscInt  k = 1;

952:     ++l;
953:     mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
954:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
955:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
956:     mpfr_set(psum, osum, MPFR_RNDN);
957:     mpfr_set(osum,  sum, MPFR_RNDN);
958:     mpfr_mul_d(h,   h,   0.5, MPFR_RNDN);
959:     mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
960:     do {
961:       mpfr_set_si(kh, k, MPFR_RNDN);
962:       mpfr_mul(kh, kh, h, MPFR_RNDN);
963:       /* Weight */
964:       mpfr_set(wk, h, MPFR_RNDN);
965:       mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
966:       mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
967:       mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
968:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
969:       mpfr_sqr(tmp, tmp, MPFR_RNDN);
970:       mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
971:       mpfr_div(wk, wk, tmp, MPFR_RNDN);
972:       /* Abscissa */
973:       mpfr_set_d(yk, 1.0, MPFR_RNDZ);
974:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
975:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
976:       mpfr_exp(tmp, msinh, MPFR_RNDN);
977:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
978:       /* Quadrature points */
979:       mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
980:       mpfr_mul(lx, lx, alpha, MPFR_RNDU);
981:       mpfr_add(lx, lx, beta, MPFR_RNDU);
982:       mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
983:       mpfr_mul(rx, rx, alpha, MPFR_RNDD);
984:       mpfr_add(rx, rx, beta, MPFR_RNDD);
985:       /* Evaluation */
986:       func(mpfr_get_d(lx, MPFR_RNDU), &lval);
987:       func(mpfr_get_d(rx, MPFR_RNDD), &rval);
988:       /* Update */
989:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
990:       mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
991:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
992:       mpfr_abs(tmp, tmp, MPFR_RNDN);
993:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
994:       mpfr_set(curTerm, tmp, MPFR_RNDN);
995:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
996:       mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
997:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
998:       mpfr_abs(tmp, tmp, MPFR_RNDN);
999:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1000:       mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
1001:       ++k;
1002:       /* Only need to evaluate every other point on refined levels */
1003:       if (l != 1) ++k;
1004:       mpfr_log10(tmp, wk, MPFR_RNDN);
1005:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1006:     } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
1007:     mpfr_sub(tmp, sum, osum, MPFR_RNDN);
1008:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1009:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1010:     d1 = mpfr_get_d(tmp, MPFR_RNDN);
1011:     mpfr_sub(tmp, sum, psum, MPFR_RNDN);
1012:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1013:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1014:     d2 = mpfr_get_d(tmp, MPFR_RNDN);
1015:     mpfr_log10(tmp, maxTerm, MPFR_RNDN);
1016:     d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
1017:     mpfr_log10(tmp, curTerm, MPFR_RNDN);
1018:     d4 = mpfr_get_d(tmp, MPFR_RNDN);
1019:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1020:   } while (d < digits && l < 8);
1021:   *sol = mpfr_get_d(sum, MPFR_RNDN);
1022:   /* Cleanup */
1023:   mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
1024:   return(0);
1025: }
1026: #else
1029: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1030: {
1031:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
1032: }
1033: #endif

1037: /* Overwrites A. Can only handle full-rank problems with m>=n
1038:  * A in column-major format
1039:  * Ainv in row-major format
1040:  * tau has length m
1041:  * worksize must be >= max(1,n)
1042:  */
1043: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1044: {
1046:   PetscBLASInt   M,N,K,lda,ldb,ldwork,info;
1047:   PetscScalar    *A,*Ainv,*R,*Q,Alpha;

1050: #if defined(PETSC_USE_COMPLEX)
1051:   {
1052:     PetscInt i,j;
1053:     PetscMalloc2(m*n,&A,m*n,&Ainv);
1054:     for (j=0; j<n; j++) {
1055:       for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1056:     }
1057:     mstride = m;
1058:   }
1059: #else
1060:   A = A_in;
1061:   Ainv = Ainv_out;
1062: #endif

1064:   PetscBLASIntCast(m,&M);
1065:   PetscBLASIntCast(n,&N);
1066:   PetscBLASIntCast(mstride,&lda);
1067:   PetscBLASIntCast(worksize,&ldwork);
1068:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1069:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1070:   PetscFPTrapPop();
1071:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1072:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1074:   /* Extract an explicit representation of Q */
1075:   Q = Ainv;
1076:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1077:   K = N;                        /* full rank */
1078:   PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1079:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

1081:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1082:   Alpha = 1.0;
1083:   ldb = lda;
1084:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1085:   /* Ainv is Q, overwritten with inverse */

1087: #if defined(PETSC_USE_COMPLEX)
1088:   {
1089:     PetscInt i;
1090:     for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1091:     PetscFree2(A,Ainv);
1092:   }
1093: #endif
1094:   return(0);
1095: }

1099: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1100: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1101: {
1103:   PetscReal      *Bv;
1104:   PetscInt       i,j;

1107:   PetscMalloc1((ninterval+1)*ndegree,&Bv);
1108:   /* Point evaluation of L_p on all the source vertices */
1109:   PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1110:   /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1111:   for (i=0; i<ninterval; i++) {
1112:     for (j=0; j<ndegree; j++) {
1113:       if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1114:       else           B[i*ndegree+j]   = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1115:     }
1116:   }
1117:   PetscFree(Bv);
1118:   return(0);
1119: }

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

1126:    Not Collective

1128:    Input Arguments:
1129: +  degree - degree of reconstruction polynomial
1130: .  nsource - number of source intervals
1131: .  sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1132: .  ntarget - number of target intervals
1133: -  targetx - sorted coordinates of target cell boundaries (length ntarget+1)

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

1138:    Level: advanced

1140: .seealso: PetscDTLegendreEval()
1141: @*/
1142: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1143: {
1145:   PetscInt       i,j,k,*bdegrees,worksize;
1146:   PetscReal      xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1147:   PetscScalar    *tau,*work;

1153:   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);
1154: #if defined(PETSC_USE_DEBUG)
1155:   for (i=0; i<nsource; i++) {
1156:     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]);
1157:   }
1158:   for (i=0; i<ntarget; i++) {
1159:     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]);
1160:   }
1161: #endif
1162:   xmin = PetscMin(sourcex[0],targetx[0]);
1163:   xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1164:   center = (xmin + xmax)/2;
1165:   hscale = (xmax - xmin)/2;
1166:   worksize = nsource;
1167:   PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1168:   PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1169:   for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1170:   for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1171:   PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1172:   PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1173:   for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1174:   PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1175:   for (i=0; i<ntarget; i++) {
1176:     PetscReal rowsum = 0;
1177:     for (j=0; j<nsource; j++) {
1178:       PetscReal sum = 0;
1179:       for (k=0; k<degree+1; k++) {
1180:         sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1181:       }
1182:       R[i*nsource+j] = sum;
1183:       rowsum += sum;
1184:     }
1185:     for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1186:   }
1187:   PetscFree4(bdegrees,sourcey,Bsource,work);
1188:   PetscFree4(tau,Bsinv,targety,Btarget);
1189:   return(0);
1190: }