Actual source code: dt.c

petsc-3.5.2 2014-09-08
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: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
 29: {

 34:   DMInitializePackage();
 35:   PetscHeaderCreate(*q,_p_PetscQuadrature,int,PETSC_OBJECT_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
 36:   (*q)->dim       = -1;
 37:   (*q)->numPoints = 0;
 38:   (*q)->points    = NULL;
 39:   (*q)->weights   = NULL;
 40:   return(0);
 41: }

 45: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
 46: {

 50:   if (!*q) return(0);
 52:   if (--((PetscObject)(*q))->refct > 0) {
 53:     *q = NULL;
 54:     return(0);
 55:   }
 56:   PetscFree((*q)->points);
 57:   PetscFree((*q)->weights);
 58:   PetscHeaderDestroy(q);
 59:   return(0);
 60: }

 64: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
 65: {
 68:   if (dim) {
 70:     *dim = q->dim;
 71:   }
 72:   if (npoints) {
 74:     *npoints = q->numPoints;
 75:   }
 76:   if (points) {
 78:     *points = q->points;
 79:   }
 80:   if (weights) {
 82:     *weights = q->weights;
 83:   }
 84:   return(0);
 85: }

 89: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
 90: {
 93:   if (dim >= 0)     q->dim       = dim;
 94:   if (npoints >= 0) q->numPoints = npoints;
 95:   if (points) {
 97:     q->points = points;
 98:   }
 99:   if (weights) {
101:     q->weights = weights;
102:   }
103:   return(0);
104: }

108: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
109: {
110:   PetscInt       q, d;

114:   PetscObjectPrintClassNamePrefixType((PetscObject)quad,viewer);
115:   PetscViewerASCIIPrintf(viewer, "Quadrature on %d points\n  (", quad->numPoints);
116:   for (q = 0; q < quad->numPoints; ++q) {
117:     for (d = 0; d < quad->dim; ++d) {
118:       if (d) PetscViewerASCIIPrintf(viewer, ", ");
119:       PetscViewerASCIIPrintf(viewer, "%g\n", (double)quad->points[q*quad->dim+d]);
120:     }
121:     PetscViewerASCIIPrintf(viewer, ") %g\n", (double)quad->weights[q]);
122:   }
123:   return(0);
124: }

128: /*@
129:    PetscDTLegendreEval - evaluate Legendre polynomial at points

131:    Not Collective

133:    Input Arguments:
134: +  npoints - number of spatial points to evaluate at
135: .  points - array of locations to evaluate at
136: .  ndegree - number of basis degrees to evaluate
137: -  degrees - sorted array of degrees to evaluate

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

144:    Level: intermediate

146: .seealso: PetscDTGaussQuadrature()
147: @*/
148: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
149: {
150:   PetscInt i,maxdegree;

153:   if (!npoints || !ndegree) return(0);
154:   maxdegree = degrees[ndegree-1];
155:   for (i=0; i<npoints; i++) {
156:     PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
157:     PetscInt  j,k;
158:     x    = points[i];
159:     pm2  = 0;
160:     pm1  = 1;
161:     pd2  = 0;
162:     pd1  = 0;
163:     pdd2 = 0;
164:     pdd1 = 0;
165:     k    = 0;
166:     if (degrees[k] == 0) {
167:       if (B) B[i*ndegree+k] = pm1;
168:       if (D) D[i*ndegree+k] = pd1;
169:       if (D2) D2[i*ndegree+k] = pdd1;
170:       k++;
171:     }
172:     for (j=1; j<=maxdegree; j++,k++) {
173:       PetscReal p,d,dd;
174:       p    = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
175:       d    = pd2 + (2*j-1)*pm1;
176:       dd   = pdd2 + (2*j-1)*pd1;
177:       pm2  = pm1;
178:       pm1  = p;
179:       pd2  = pd1;
180:       pd1  = d;
181:       pdd2 = pdd1;
182:       pdd1 = dd;
183:       if (degrees[k] == j) {
184:         if (B) B[i*ndegree+k] = p;
185:         if (D) D[i*ndegree+k] = d;
186:         if (D2) D2[i*ndegree+k] = dd;
187:       }
188:     }
189:   }
190:   return(0);
191: }

195: /*@
196:    PetscDTGaussQuadrature - create Gauss quadrature

198:    Not Collective

200:    Input Arguments:
201: +  npoints - number of points
202: .  a - left end of interval (often-1)
203: -  b - right end of interval (often +1)

205:    Output Arguments:
206: +  x - quadrature points
207: -  w - quadrature weights

209:    Level: intermediate

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

214: .seealso: PetscDTLegendreEval()
215: @*/
216: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
217: {
219:   PetscInt       i;
220:   PetscReal      *work;
221:   PetscScalar    *Z;
222:   PetscBLASInt   N,LDZ,info;

225:   PetscCitationsRegister(GaussCitation, &GaussCite);
226:   /* Set up the Golub-Welsch system */
227:   for (i=0; i<npoints; i++) {
228:     x[i] = 0;                   /* diagonal is 0 */
229:     if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
230:   }
231:   PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
232:   PetscBLASIntCast(npoints,&N);
233:   LDZ  = N;
234:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
235:   PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
236:   PetscFPTrapPop();
237:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");

239:   for (i=0; i<(npoints+1)/2; i++) {
240:     PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
241:     x[i]           = (a+b)/2 - y*(b-a)/2;
242:     x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;

244:     w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
245:   }
246:   PetscFree2(Z,work);
247:   return(0);
248: }

252: /*@
253:   PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature

255:   Not Collective

257:   Input Arguments:
258: + dim     - The spatial dimension
259: . npoints - number of points in one dimension
260: . a       - left end of interval (often-1)
261: - b       - right end of interval (often +1)

263:   Output Argument:
264: . q - A PetscQuadrature object

266:   Level: intermediate

268: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
269: @*/
270: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
271: {
272:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k;
273:   PetscReal     *x, *w, *xw, *ww;

277:   PetscMalloc1(totpoints*dim,&x);
278:   PetscMalloc1(totpoints,&w);
279:   /* Set up the Golub-Welsch system */
280:   switch (dim) {
281:   case 0:
282:     PetscFree(x);
283:     PetscFree(w);
284:     PetscMalloc1(1, &x);
285:     PetscMalloc1(1, &w);
286:     x[0] = 0.0;
287:     w[0] = 1.0;
288:     break;
289:   case 1:
290:     PetscDTGaussQuadrature(npoints, a, b, x, w);
291:     break;
292:   case 2:
293:     PetscMalloc2(npoints,&xw,npoints,&ww);
294:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
295:     for (i = 0; i < npoints; ++i) {
296:       for (j = 0; j < npoints; ++j) {
297:         x[(i*npoints+j)*dim+0] = xw[i];
298:         x[(i*npoints+j)*dim+1] = xw[j];
299:         w[i*npoints+j]         = ww[i] * ww[j];
300:       }
301:     }
302:     PetscFree2(xw,ww);
303:     break;
304:   case 3:
305:     PetscMalloc2(npoints,&xw,npoints,&ww);
306:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
307:     for (i = 0; i < npoints; ++i) {
308:       for (j = 0; j < npoints; ++j) {
309:         for (k = 0; k < npoints; ++k) {
310:           x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
311:           x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
312:           x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
313:           w[(i*npoints+j)*npoints+k]         = ww[i] * ww[j] * ww[k];
314:         }
315:       }
316:     }
317:     PetscFree2(xw,ww);
318:     break;
319:   default:
320:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
321:   }
322:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
323:   PetscQuadratureSetData(*q, dim, totpoints, x, w);
324:   return(0);
325: }

329: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
330:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
331: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
332: {
333:   PetscReal f = 1.0;
334:   PetscInt  i;

337:   for (i = 1; i < n+1; ++i) f *= i;
338:   *factorial = f;
339:   return(0);
340: }

344: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
345:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
346: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
347: {
348:   PetscReal apb, pn1, pn2;
349:   PetscInt  k;

352:   if (!n) {*P = 1.0; return(0);}
353:   if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
354:   apb = a + b;
355:   pn2 = 1.0;
356:   pn1 = 0.5 * (a - b + (apb + 2.0) * x);
357:   *P  = 0.0;
358:   for (k = 2; k < n+1; ++k) {
359:     PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
360:     PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
361:     PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
362:     PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);

364:     a2  = a2 / a1;
365:     a3  = a3 / a1;
366:     a4  = a4 / a1;
367:     *P  = (a2 + a3 * x) * pn1 - a4 * pn2;
368:     pn2 = pn1;
369:     pn1 = *P;
370:   }
371:   return(0);
372: }

376: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
377: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
378: {
379:   PetscReal      nP;

383:   if (!n) {*P = 0.0; return(0);}
384:   PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
385:   *P   = 0.5 * (a + b + n + 1) * nP;
386:   return(0);
387: }

391: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
392: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
393: {
395:   *xi  = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
396:   *eta = y;
397:   return(0);
398: }

402: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
403: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
404: {
406:   *xi   = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
407:   *eta  = 0.5  * (1.0 + y) * (1.0 - z) - 1.0;
408:   *zeta = z;
409:   return(0);
410: }

414: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
415: {
416:   PetscInt       maxIter = 100;
417:   PetscReal      eps     = 1.0e-8;
418:   PetscReal      a1, a2, a3, a4, a5, a6;
419:   PetscInt       k;


424:   a1      = PetscPowReal(2.0, a+b+1);
425: #if defined(PETSC_HAVE_TGAMMA)
426:   a2      = PetscTGamma(a + npoints + 1);
427:   a3      = PetscTGamma(b + npoints + 1);
428:   a4      = PetscTGamma(a + b + npoints + 1);
429: #else
430:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
431: #endif

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

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

446:       for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
447:       PetscDTComputeJacobi(a, b, npoints, r, &f);
448:       PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
449:       delta = f / (fp - f * s);
450:       r     = r - delta;
451:       if (PetscAbsReal(delta) < eps) break;
452:     }
453:     x[k] = r;
454:     PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
455:     w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
456:   }
457:   return(0);
458: }

462: /*@C
463:   PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex

465:   Not Collective

467:   Input Arguments:
468: + dim   - The simplex dimension
469: . order - The number of points in one dimension
470: . a     - left end of interval (often-1)
471: - b     - right end of interval (often +1)

473:   Output Argument:
474: . q - A PetscQuadrature object

476:   Level: intermediate

478:   References:
479:   Karniadakis and Sherwin.
480:   FIAT

482: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
483: @*/
484: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt order, PetscReal a, PetscReal b, PetscQuadrature *q)
485: {
486:   PetscInt       npoints = dim > 1 ? dim > 2 ? order*PetscSqr(order) : PetscSqr(order) : order;
487:   PetscReal     *px, *wx, *py, *wy, *pz, *wz, *x, *w;
488:   PetscInt       i, j, k;

492:   if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
493:   PetscMalloc1(npoints*dim, &x);
494:   PetscMalloc1(npoints, &w);
495:   switch (dim) {
496:   case 0:
497:     PetscFree(x);
498:     PetscFree(w);
499:     PetscMalloc1(1, &x);
500:     PetscMalloc1(1, &w);
501:     x[0] = 0.0;
502:     w[0] = 1.0;
503:     break;
504:   case 1:
505:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, x, w);
506:     break;
507:   case 2:
508:     PetscMalloc4(order,&px,order,&wx,order,&py,order,&wy);
509:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
510:     PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
511:     for (i = 0; i < order; ++i) {
512:       for (j = 0; j < order; ++j) {
513:         PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*order+j)*2+0], &x[(i*order+j)*2+1]);
514:         w[i*order+j] = 0.5 * wx[i] * wy[j];
515:       }
516:     }
517:     PetscFree4(px,wx,py,wy);
518:     break;
519:   case 3:
520:     PetscMalloc6(order,&px,order,&wx,order,&py,order,&wy,order,&pz,order,&wz);
521:     PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);
522:     PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);
523:     PetscDTGaussJacobiQuadrature1D_Internal(order, 2.0, 0.0, pz, wz);
524:     for (i = 0; i < order; ++i) {
525:       for (j = 0; j < order; ++j) {
526:         for (k = 0; k < order; ++k) {
527:           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]);
528:           w[(i*order+j)*order+k] = 0.125 * wx[i] * wy[j] * wz[k];
529:         }
530:       }
531:     }
532:     PetscFree6(px,wx,py,wy,pz,wz);
533:     break;
534:   default:
535:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
536:   }
537:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
538:   PetscQuadratureSetData(*q, dim, npoints, x, w);
539:   return(0);
540: }

544: /* Overwrites A. Can only handle full-rank problems with m>=n
545:  * A in column-major format
546:  * Ainv in row-major format
547:  * tau has length m
548:  * worksize must be >= max(1,n)
549:  */
550: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
551: {
553:   PetscBLASInt   M,N,K,lda,ldb,ldwork,info;
554:   PetscScalar    *A,*Ainv,*R,*Q,Alpha;

557: #if defined(PETSC_USE_COMPLEX)
558:   {
559:     PetscInt i,j;
560:     PetscMalloc2(m*n,&A,m*n,&Ainv);
561:     for (j=0; j<n; j++) {
562:       for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
563:     }
564:     mstride = m;
565:   }
566: #else
567:   A = A_in;
568:   Ainv = Ainv_out;
569: #endif

571:   PetscBLASIntCast(m,&M);
572:   PetscBLASIntCast(n,&N);
573:   PetscBLASIntCast(mstride,&lda);
574:   PetscBLASIntCast(worksize,&ldwork);
575:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
576:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
577:   PetscFPTrapPop();
578:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
579:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

581:   /* Extract an explicit representation of Q */
582:   Q = Ainv;
583:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
584:   K = N;                        /* full rank */
585:   PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
586:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

588:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
589:   Alpha = 1.0;
590:   ldb = lda;
591:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
592:   /* Ainv is Q, overwritten with inverse */

594: #if defined(PETSC_USE_COMPLEX)
595:   {
596:     PetscInt i;
597:     for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
598:     PetscFree2(A,Ainv);
599:   }
600: #endif
601:   return(0);
602: }

606: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
607: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
608: {
610:   PetscReal      *Bv;
611:   PetscInt       i,j;

614:   PetscMalloc1((ninterval+1)*ndegree,&Bv);
615:   /* Point evaluation of L_p on all the source vertices */
616:   PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
617:   /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
618:   for (i=0; i<ninterval; i++) {
619:     for (j=0; j<ndegree; j++) {
620:       if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
621:       else           B[i*ndegree+j]   = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
622:     }
623:   }
624:   PetscFree(Bv);
625:   return(0);
626: }

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

633:    Not Collective

635:    Input Arguments:
636: +  degree - degree of reconstruction polynomial
637: .  nsource - number of source intervals
638: .  sourcex - sorted coordinates of source cell boundaries (length nsource+1)
639: .  ntarget - number of target intervals
640: -  targetx - sorted coordinates of target cell boundaries (length ntarget+1)

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

645:    Level: advanced

647: .seealso: PetscDTLegendreEval()
648: @*/
649: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
650: {
652:   PetscInt       i,j,k,*bdegrees,worksize;
653:   PetscReal      xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
654:   PetscScalar    *tau,*work;

660:   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);
661: #if defined(PETSC_USE_DEBUG)
662:   for (i=0; i<nsource; i++) {
663:     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]);
664:   }
665:   for (i=0; i<ntarget; i++) {
666:     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]);
667:   }
668: #endif
669:   xmin = PetscMin(sourcex[0],targetx[0]);
670:   xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
671:   center = (xmin + xmax)/2;
672:   hscale = (xmax - xmin)/2;
673:   worksize = nsource;
674:   PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
675:   PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
676:   for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
677:   for (i=0; i<=degree; i++) bdegrees[i] = i+1;
678:   PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
679:   PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
680:   for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
681:   PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
682:   for (i=0; i<ntarget; i++) {
683:     PetscReal rowsum = 0;
684:     for (j=0; j<nsource; j++) {
685:       PetscReal sum = 0;
686:       for (k=0; k<degree+1; k++) {
687:         sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
688:       }
689:       R[i*nsource+j] = sum;
690:       rowsum += sum;
691:     }
692:     for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
693:   }
694:   PetscFree4(bdegrees,sourcey,Bsource,work);
695:   PetscFree4(tau,Bsinv,targety,Btarget);
696:   return(0);
697: }