Actual source code: dt.c

petsc-master 2015-01-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

  8: #include <petscdt.h>            /*I "petscdt.h" I*/
  9: #include <petscblaslapack.h>
 10: #include <petsc-private/petscimpl.h>
 11: #include <petsc-private/dtimpl.h>
 12: #include <petscviewer.h>
 13: #include <petscdmplex.h>
 14: #include <petscdmshell.h>

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

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

 31:   Collective on MPI_Comm

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

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

 39:   Level: beginner

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

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

 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, 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, &Nq, &points, &weights);
 91:   PetscMalloc1(Nq*dim, &p);
 92:   PetscMalloc1(Nq, &w);
 93:   PetscMemcpy(p, points, Nq*dim * sizeof(PetscReal));
 94:   PetscMemcpy(w, weights, Nq * sizeof(PetscReal));
 95:   PetscQuadratureSetData(*r, dim, Nq, p, w);
 96:   return(0);
 97: }

101: /*@
102:   PetscQuadratureDestroy - Destroys a PetscQuadrature object

104:   Collective on PetscQuadrature

106:   Input Parameter:
107: . q  - The PetscQuadrature object

109:   Level: beginner

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

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

133: /*@
134:   PetscQuadratureGetOrder - Return the quadrature information

136:   Not collective

138:   Input Parameter:
139: . q - The PetscQuadrature object

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

144:   Output Parameter:

146:   Level: intermediate

148: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
149: @*/
150: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
151: {
155:   *order = q->order;
156:   return(0);
157: }

161: /*@
162:   PetscQuadratureSetOrder - Return the quadrature information

164:   Not collective

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

170:   Level: intermediate

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

184: /*@C
185:   PetscQuadratureGetData - Returns the data defining the quadrature

187:   Not collective

189:   Input Parameter:
190: . q  - The PetscQuadrature object

192:   Output Parameters:
193: + dim - The spatial dimension
194: . npoints - The number of quadrature points
195: . points - The coordinates of each quadrature point
196: - weights - The weight of each quadrature point

198:   Level: intermediate

200: .keywords: PetscQuadrature, quadrature
201: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
202: @*/
203: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
204: {
207:   if (dim) {
209:     *dim = q->dim;
210:   }
211:   if (npoints) {
213:     *npoints = q->numPoints;
214:   }
215:   if (points) {
217:     *points = q->points;
218:   }
219:   if (weights) {
221:     *weights = q->weights;
222:   }
223:   return(0);
224: }

228: /*@C
229:   PetscQuadratureSetData - Sets the data defining the quadrature

231:   Not collective

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

240:   Level: intermediate

242: .keywords: PetscQuadrature, quadrature
243: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
244: @*/
245: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
246: {
249:   if (dim >= 0)     q->dim       = dim;
250:   if (npoints >= 0) q->numPoints = npoints;
251:   if (points) {
253:     q->points = points;
254:   }
255:   if (weights) {
257:     q->weights = weights;
258:   }
259:   return(0);
260: }

264: /*@C
265:   PetscQuadratureView - Views a PetscQuadrature object

267:   Collective on PetscQuadrature

269:   Input Parameters:
270: + q  - The PetscQuadrature object
271: - viewer - The PetscViewer object

273:   Level: beginner

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

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

298: /*@
299:    PetscDTLegendreEval - evaluate Legendre polynomial at points

301:    Not Collective

303:    Input Arguments:
304: +  npoints - number of spatial points to evaluate at
305: .  points - array of locations to evaluate at
306: .  ndegree - number of basis degrees to evaluate
307: -  degrees - sorted array of degrees to evaluate

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

314:    Level: intermediate

316: .seealso: PetscDTGaussQuadrature()
317: @*/
318: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
319: {
320:   PetscInt i,maxdegree;

323:   if (!npoints || !ndegree) return(0);
324:   maxdegree = degrees[ndegree-1];
325:   for (i=0; i<npoints; i++) {
326:     PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
327:     PetscInt  j,k;
328:     x    = points[i];
329:     pm2  = 0;
330:     pm1  = 1;
331:     pd2  = 0;
332:     pd1  = 0;
333:     pdd2 = 0;
334:     pdd1 = 0;
335:     k    = 0;
336:     if (degrees[k] == 0) {
337:       if (B) B[i*ndegree+k] = pm1;
338:       if (D) D[i*ndegree+k] = pd1;
339:       if (D2) D2[i*ndegree+k] = pdd1;
340:       k++;
341:     }
342:     for (j=1; j<=maxdegree; j++,k++) {
343:       PetscReal p,d,dd;
344:       p    = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
345:       d    = pd2 + (2*j-1)*pm1;
346:       dd   = pdd2 + (2*j-1)*pd1;
347:       pm2  = pm1;
348:       pm1  = p;
349:       pd2  = pd1;
350:       pd1  = d;
351:       pdd2 = pdd1;
352:       pdd1 = dd;
353:       if (degrees[k] == j) {
354:         if (B) B[i*ndegree+k] = p;
355:         if (D) D[i*ndegree+k] = d;
356:         if (D2) D2[i*ndegree+k] = dd;
357:       }
358:     }
359:   }
360:   return(0);
361: }

365: /*@
366:    PetscDTGaussQuadrature - create Gauss quadrature

368:    Not Collective

370:    Input Arguments:
371: +  npoints - number of points
372: .  a - left end of interval (often-1)
373: -  b - right end of interval (often +1)

375:    Output Arguments:
376: +  x - quadrature points
377: -  w - quadrature weights

379:    Level: intermediate

381:    References:
382:    Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 221--230, 1969.

384: .seealso: PetscDTLegendreEval()
385: @*/
386: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
387: {
389:   PetscInt       i;
390:   PetscReal      *work;
391:   PetscScalar    *Z;
392:   PetscBLASInt   N,LDZ,info;

395:   PetscCitationsRegister(GaussCitation, &GaussCite);
396:   /* Set up the Golub-Welsch system */
397:   for (i=0; i<npoints; i++) {
398:     x[i] = 0;                   /* diagonal is 0 */
399:     if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
400:   }
401:   PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
402:   PetscBLASIntCast(npoints,&N);
403:   LDZ  = N;
404:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
405:   PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
406:   PetscFPTrapPop();
407:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");

409:   for (i=0; i<(npoints+1)/2; i++) {
410:     PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
411:     x[i]           = (a+b)/2 - y*(b-a)/2;
412:     x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;

414:     w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
415:   }
416:   PetscFree2(Z,work);
417:   return(0);
418: }

422: /*@
423:   PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature

425:   Not Collective

427:   Input Arguments:
428: + dim     - The spatial dimension
429: . npoints - number of points in one dimension
430: . a       - left end of interval (often-1)
431: - b       - right end of interval (often +1)

433:   Output Argument:
434: . q - A PetscQuadrature object

436:   Level: intermediate

438: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
439: @*/
440: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
441: {
442:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k;
443:   PetscReal     *x, *w, *xw, *ww;

447:   PetscMalloc1(totpoints*dim,&x);
448:   PetscMalloc1(totpoints,&w);
449:   /* Set up the Golub-Welsch system */
450:   switch (dim) {
451:   case 0:
452:     PetscFree(x);
453:     PetscFree(w);
454:     PetscMalloc1(1, &x);
455:     PetscMalloc1(1, &w);
456:     x[0] = 0.0;
457:     w[0] = 1.0;
458:     break;
459:   case 1:
460:     PetscDTGaussQuadrature(npoints, a, b, x, w);
461:     break;
462:   case 2:
463:     PetscMalloc2(npoints,&xw,npoints,&ww);
464:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
465:     for (i = 0; i < npoints; ++i) {
466:       for (j = 0; j < npoints; ++j) {
467:         x[(i*npoints+j)*dim+0] = xw[i];
468:         x[(i*npoints+j)*dim+1] = xw[j];
469:         w[i*npoints+j]         = ww[i] * ww[j];
470:       }
471:     }
472:     PetscFree2(xw,ww);
473:     break;
474:   case 3:
475:     PetscMalloc2(npoints,&xw,npoints,&ww);
476:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
477:     for (i = 0; i < npoints; ++i) {
478:       for (j = 0; j < npoints; ++j) {
479:         for (k = 0; k < npoints; ++k) {
480:           x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
481:           x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
482:           x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
483:           w[(i*npoints+j)*npoints+k]         = ww[i] * ww[j] * ww[k];
484:         }
485:       }
486:     }
487:     PetscFree2(xw,ww);
488:     break;
489:   default:
490:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
491:   }
492:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
493:   PetscQuadratureSetOrder(*q, npoints);
494:   PetscQuadratureSetData(*q, dim, totpoints, x, w);
495:   return(0);
496: }

500: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
501:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
502: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
503: {
504:   PetscReal f = 1.0;
505:   PetscInt  i;

508:   for (i = 1; i < n+1; ++i) f *= i;
509:   *factorial = f;
510:   return(0);
511: }

515: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
516:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
517: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
518: {
519:   PetscReal apb, pn1, pn2;
520:   PetscInt  k;

523:   if (!n) {*P = 1.0; return(0);}
524:   if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
525:   apb = a + b;
526:   pn2 = 1.0;
527:   pn1 = 0.5 * (a - b + (apb + 2.0) * x);
528:   *P  = 0.0;
529:   for (k = 2; k < n+1; ++k) {
530:     PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
531:     PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
532:     PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
533:     PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);

535:     a2  = a2 / a1;
536:     a3  = a3 / a1;
537:     a4  = a4 / a1;
538:     *P  = (a2 + a3 * x) * pn1 - a4 * pn2;
539:     pn2 = pn1;
540:     pn1 = *P;
541:   }
542:   return(0);
543: }

547: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
548: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
549: {
550:   PetscReal      nP;

554:   if (!n) {*P = 0.0; return(0);}
555:   PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
556:   *P   = 0.5 * (a + b + n + 1) * nP;
557:   return(0);
558: }

562: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
563: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
564: {
566:   *xi  = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
567:   *eta = y;
568:   return(0);
569: }

573: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
574: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
575: {
577:   *xi   = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
578:   *eta  = 0.5  * (1.0 + y) * (1.0 - z) - 1.0;
579:   *zeta = z;
580:   return(0);
581: }

585: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
586: {
587:   PetscInt       maxIter = 100;
588:   PetscReal      eps     = 1.0e-8;
589:   PetscReal      a1, a2, a3, a4, a5, a6;
590:   PetscInt       k;


595:   a1      = PetscPowReal(2.0, a+b+1);
596: #if defined(PETSC_HAVE_TGAMMA)
597:   a2      = PetscTGamma(a + npoints + 1);
598:   a3      = PetscTGamma(b + npoints + 1);
599:   a4      = PetscTGamma(a + b + npoints + 1);
600: #else
601:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
602: #endif

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

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

617:       for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
618:       PetscDTComputeJacobi(a, b, npoints, r, &f);
619:       PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
620:       delta = f / (fp - f * s);
621:       r     = r - delta;
622:       if (PetscAbsReal(delta) < eps) break;
623:     }
624:     x[k] = r;
625:     PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
626:     w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
627:   }
628:   return(0);
629: }

633: /*@C
634:   PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex

636:   Not Collective

638:   Input Arguments:
639: + dim   - The simplex dimension
640: . order - The number of points in one dimension
641: . a     - left end of interval (often-1)
642: - b     - right end of interval (often +1)

644:   Output Argument:
645: . q - A PetscQuadrature object

647:   Level: intermediate

649:   References:
650:   Karniadakis and Sherwin.
651:   FIAT

653: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
654: @*/
655: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt order, PetscReal a, PetscReal b, PetscQuadrature *q)
656: {
657:   PetscInt       npoints = dim > 1 ? dim > 2 ? order*PetscSqr(order) : PetscSqr(order) : order;
658:   PetscReal     *px, *wx, *py, *wy, *pz, *wz, *x, *w;
659:   PetscInt       i, j, k;

663:   if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
664:   PetscMalloc1(npoints*dim, &x);
665:   PetscMalloc1(npoints, &w);
666:   switch (dim) {
667:   case 0:
668:     PetscFree(x);
669:     PetscFree(w);
670:     PetscMalloc1(1, &x);
671:     PetscMalloc1(1, &w);
672:     x[0] = 0.0;
673:     w[0] = 1.0;
674:     break;
675:   case 1:
676:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, x, w);
677:     break;
678:   case 2:
679:     PetscMalloc4(order,&px,order,&wx,order,&py,order,&wy);
680:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
681:     PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
682:     for (i = 0; i < order; ++i) {
683:       for (j = 0; j < order; ++j) {
684:         PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*order+j)*2+0], &x[(i*order+j)*2+1]);
685:         w[i*order+j] = 0.5 * wx[i] * wy[j];
686:       }
687:     }
688:     PetscFree4(px,wx,py,wy);
689:     break;
690:   case 3:
691:     PetscMalloc6(order,&px,order,&wx,order,&py,order,&wy,order,&pz,order,&wz);
692:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
693:     PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
694:     PetscDTGaussJacobiQuadrature1D_Internal(order, 2.0, 0.0, pz, wz);
695:     for (i = 0; i < order; ++i) {
696:       for (j = 0; j < order; ++j) {
697:         for (k = 0; k < order; ++k) {
698:           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]);
699:           w[(i*order+j)*order+k] = 0.125 * wx[i] * wy[j] * wz[k];
700:         }
701:       }
702:     }
703:     PetscFree6(px,wx,py,wy,pz,wz);
704:     break;
705:   default:
706:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
707:   }
708:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
709:   PetscQuadratureSetOrder(*q, order);
710:   PetscQuadratureSetData(*q, dim, npoints, x, w);
711:   return(0);
712: }

716: /* Overwrites A. Can only handle full-rank problems with m>=n
717:  * A in column-major format
718:  * Ainv in row-major format
719:  * tau has length m
720:  * worksize must be >= max(1,n)
721:  */
722: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
723: {
725:   PetscBLASInt   M,N,K,lda,ldb,ldwork,info;
726:   PetscScalar    *A,*Ainv,*R,*Q,Alpha;

729: #if defined(PETSC_USE_COMPLEX)
730:   {
731:     PetscInt i,j;
732:     PetscMalloc2(m*n,&A,m*n,&Ainv);
733:     for (j=0; j<n; j++) {
734:       for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
735:     }
736:     mstride = m;
737:   }
738: #else
739:   A = A_in;
740:   Ainv = Ainv_out;
741: #endif

743:   PetscBLASIntCast(m,&M);
744:   PetscBLASIntCast(n,&N);
745:   PetscBLASIntCast(mstride,&lda);
746:   PetscBLASIntCast(worksize,&ldwork);
747:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
748:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
749:   PetscFPTrapPop();
750:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
751:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

753:   /* Extract an explicit representation of Q */
754:   Q = Ainv;
755:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
756:   K = N;                        /* full rank */
757:   PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
758:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

760:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
761:   Alpha = 1.0;
762:   ldb = lda;
763:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
764:   /* Ainv is Q, overwritten with inverse */

766: #if defined(PETSC_USE_COMPLEX)
767:   {
768:     PetscInt i;
769:     for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
770:     PetscFree2(A,Ainv);
771:   }
772: #endif
773:   return(0);
774: }

778: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
779: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
780: {
782:   PetscReal      *Bv;
783:   PetscInt       i,j;

786:   PetscMalloc1((ninterval+1)*ndegree,&Bv);
787:   /* Point evaluation of L_p on all the source vertices */
788:   PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
789:   /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
790:   for (i=0; i<ninterval; i++) {
791:     for (j=0; j<ndegree; j++) {
792:       if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
793:       else           B[i*ndegree+j]   = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
794:     }
795:   }
796:   PetscFree(Bv);
797:   return(0);
798: }

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

805:    Not Collective

807:    Input Arguments:
808: +  degree - degree of reconstruction polynomial
809: .  nsource - number of source intervals
810: .  sourcex - sorted coordinates of source cell boundaries (length nsource+1)
811: .  ntarget - number of target intervals
812: -  targetx - sorted coordinates of target cell boundaries (length ntarget+1)

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

817:    Level: advanced

819: .seealso: PetscDTLegendreEval()
820: @*/
821: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
822: {
824:   PetscInt       i,j,k,*bdegrees,worksize;
825:   PetscReal      xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
826:   PetscScalar    *tau,*work;

832:   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);
833: #if defined(PETSC_USE_DEBUG)
834:   for (i=0; i<nsource; i++) {
835:     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]);
836:   }
837:   for (i=0; i<ntarget; i++) {
838:     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]);
839:   }
840: #endif
841:   xmin = PetscMin(sourcex[0],targetx[0]);
842:   xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
843:   center = (xmin + xmax)/2;
844:   hscale = (xmax - xmin)/2;
845:   worksize = nsource;
846:   PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
847:   PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
848:   for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
849:   for (i=0; i<=degree; i++) bdegrees[i] = i+1;
850:   PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
851:   PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
852:   for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
853:   PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
854:   for (i=0; i<ntarget; i++) {
855:     PetscReal rowsum = 0;
856:     for (j=0; j<nsource; j++) {
857:       PetscReal sum = 0;
858:       for (k=0; k<degree+1; k++) {
859:         sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
860:       }
861:       R[i*nsource+j] = sum;
862:       rowsum += sum;
863:     }
864:     for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
865:   }
866:   PetscFree4(bdegrees,sourcey,Bsource,work);
867:   PetscFree4(tau,Bsinv,targety,Btarget);
868:   return(0);
869: }