Actual source code: dt.c

petsc-master 2019-12-13
Report Typos and Errors
  1: /* Discretization tools */

  3:  #include <petscdt.h>
  4:  #include <petscblaslapack.h>
  5:  #include <petsc/private/petscimpl.h>
  6:  #include <petsc/private/dtimpl.h>
  7:  #include <petscviewer.h>
  8:  #include <petscdmplex.h>
  9:  #include <petscdmshell.h>

 11: #if defined(PETSC_HAVE_MPFR)
 12: #include <mpfr.h>
 13: #endif

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


 26: PetscClassId PETSCQUADRATURE_CLASSID = 0;

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

 31:   Collective

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

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

 39:   Level: beginner

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

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

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

 63:   Collective on q

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

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

 71:   Level: beginner

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

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

 96: /*@
 97:   PetscQuadratureDestroy - Destroys a PetscQuadrature object

 99:   Collective on q

101:   Input Parameter:
102: . q  - The PetscQuadrature object

104:   Level: beginner

106: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
107: @*/
108: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
109: {

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

125: /*@
126:   PetscQuadratureGetOrder - Return the order of the method

128:   Not collective

130:   Input Parameter:
131: . q - The PetscQuadrature object

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

136:   Level: intermediate

138: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
139: @*/
140: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
141: {
145:   *order = q->order;
146:   return(0);
147: }

149: /*@
150:   PetscQuadratureSetOrder - Return the order of the method

152:   Not collective

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

158:   Level: intermediate

160: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
161: @*/
162: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
163: {
166:   q->order = order;
167:   return(0);
168: }

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

173:   Not collective

175:   Input Parameter:
176: . q - The PetscQuadrature object

178:   Output Parameter:
179: . Nc - The number of components

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

183:   Level: intermediate

185: .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
186: @*/
187: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
188: {
192:   *Nc = q->Nc;
193:   return(0);
194: }

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

199:   Not collective

201:   Input Parameters:
202: + q  - The PetscQuadrature object
203: - Nc - The number of components

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

207:   Level: intermediate

209: .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
210: @*/
211: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
212: {
215:   q->Nc = Nc;
216:   return(0);
217: }

219: /*@C
220:   PetscQuadratureGetData - Returns the data defining the quadrature

222:   Not collective

224:   Input Parameter:
225: . q  - The PetscQuadrature object

227:   Output Parameters:
228: + dim - The spatial dimension
229: . Nc - The number of components
230: . npoints - The number of quadrature points
231: . points - The coordinates of each quadrature point
232: - weights - The weight of each quadrature point

234:   Level: intermediate

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

239: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
240: @*/
241: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
242: {
245:   if (dim) {
247:     *dim = q->dim;
248:   }
249:   if (Nc) {
251:     *Nc = q->Nc;
252:   }
253:   if (npoints) {
255:     *npoints = q->numPoints;
256:   }
257:   if (points) {
259:     *points = q->points;
260:   }
261:   if (weights) {
263:     *weights = q->weights;
264:   }
265:   return(0);
266: }

268: static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[])
269: {
270:   PetscScalar    *Js, *Jinvs;
271:   PetscInt       i, j, k;
272:   PetscBLASInt   bm, bn, info;

276:   PetscBLASIntCast(m, &bm);
277:   PetscBLASIntCast(n, &bn);
278: #if defined(PETSC_USE_COMPLEX)
279:   PetscMalloc2(m*n, &Js, m*n, &Jinvs);
280:   for (i = 0; i < m*n; i++) Js[i] = J[i];
281: #else
282:   Js = (PetscReal *) J;
283:   Jinvs = Jinv;
284: #endif
285:   if (m == n) {
286:     PetscBLASInt *pivots;
287:     PetscScalar *W;

289:     PetscMalloc2(m, &pivots, m, &W);

291:     PetscArraycpy(Jinvs, Js, m * m);
292:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info));
293:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
294:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &info));
295:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
296:     PetscFree2(pivots, W);
297:   } else if (m < n) {
298:     PetscScalar *JJT;
299:     PetscBLASInt *pivots;
300:     PetscScalar *W;

302:     PetscMalloc1(m*m, &JJT);
303:     PetscMalloc2(m, &pivots, m, &W);
304:     for (i = 0; i < m; i++) {
305:       for (j = 0; j < m; j++) {
306:         PetscScalar val = 0.;

308:         for (k = 0; k < n; k++) val += Js[i * n + k] * Js[j * n + k];
309:         JJT[i * m + j] = val;
310:       }
311:     }

313:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info));
314:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
315:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info));
316:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
317:     for (i = 0; i < n; i++) {
318:       for (j = 0; j < m; j++) {
319:         PetscScalar val = 0.;

321:         for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j];
322:         Jinvs[i * m + j] = val;
323:       }
324:     }
325:     PetscFree2(pivots, W);
326:     PetscFree(JJT);
327:   } else {
328:     PetscScalar *JTJ;
329:     PetscBLASInt *pivots;
330:     PetscScalar *W;

332:     PetscMalloc1(n*n, &JTJ);
333:     PetscMalloc2(n, &pivots, n, &W);
334:     for (i = 0; i < n; i++) {
335:       for (j = 0; j < n; j++) {
336:         PetscScalar val = 0.;

338:         for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j];
339:         JTJ[i * n + j] = val;
340:       }
341:     }

343:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bm, pivots, &info));
344:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
345:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info));
346:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
347:     for (i = 0; i < n; i++) {
348:       for (j = 0; j < m; j++) {
349:         PetscScalar val = 0.;

351:         for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k];
352:         Jinvs[i * m + j] = val;
353:       }
354:     }
355:     PetscFree2(pivots, W);
356:     PetscFree(JTJ);
357:   }
358: #if defined(PETSC_USE_COMPLEX)
359:   for (i = 0; i < m*n; i++) Jinv[i] = PetscRealPart(Jinvs[i]);
360:   PetscFree2(Js, Jinvs);
361: #endif
362:   return(0);
363: }

365: /*@
366:    PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation.

368:    Collecive on PetscQuadrature

370:    Input Arguments:
371: +  q - the quadrature functional
372: .  imageDim - the dimension of the image of the transformation
373: .  origin - a point in the original space
374: .  originImage - the image of the origin under the transformation
375: .  J - the Jacobian of the image: an [imageDim x dim] matrix in row major order
376: -  formDegree - transform the quadrature weights as k-forms of this form degree (if the number of components is a multiple of (dim choose formDegree), it is assumed that they represent multiple k-forms) [see PetscDTAltVPullback() for interpretation of formDegree]

378:    Output Arguments:
379: .  Jinvstarq - a quadrature rule where each point is the image of a point in the original quadrature rule, and where the k-form weights have been pulled-back by the pseudoinverse of J to the k-form weights in the image space.

381:    Note: the new quadrature rule will have a different number of components if spaces have different dimensions.  For example, pushing a 2-form forward from a two dimensional space to a three dimensional space changes the number of components from 1 to 3.

383: .seealso: PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
384: @*/
385: PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq)
386: {
387:   PetscInt         dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c;
388:   const PetscReal *points;
389:   const PetscReal *weights;
390:   PetscReal       *imagePoints, *imageWeights;
391:   PetscReal       *Jinv;
392:   PetscReal       *Jinvstar;
393:   PetscErrorCode   ierr;

397:   if (imageDim < PetscAbsInt(formDegree)) SETERRQ2(PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Cannot represent a %D-form in %D dimensions", PetscAbsInt(formDegree), imageDim);
398:   PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights);
399:   PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize);
400:   if (Nc % formSize) SETERRQ2(PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Number of components %D is not a multiple of formSize %D\n", Nc, formSize);
401:   Ncopies = Nc / formSize;
402:   PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize);
403:   imageNc = Ncopies * imageFormSize;
404:   PetscMalloc1(Npoints * imageDim, &imagePoints);
405:   PetscMalloc1(Npoints * imageNc, &imageWeights);
406:   PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar);
407:   PetscDTJacobianInverse_Internal(dim, imageDim, J, Jinv);
408:   PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar);
409:   for (pt = 0; pt < Npoints; pt++) {
410:     const PetscReal *point = &points[pt * dim];
411:     PetscReal       *imagePoint = &imagePoints[pt * imageDim];

413:     for (i = 0; i < imageDim; i++) {
414:       PetscReal val = originImage[i];

416:       for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]);
417:       imagePoint[i] = val;
418:     }
419:     for (c = 0; c < Ncopies; c++) {
420:       const PetscReal *form = &weights[pt * Nc + c * formSize];
421:       PetscReal       *imageForm = &imageWeights[pt * imageNc + c * imageFormSize];

423:       for (i = 0; i < imageFormSize; i++) {
424:         PetscReal val = 0.;

426:         for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j];
427:         imageForm[i] = val;
428:       }
429:     }
430:   }
431:   PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq);
432:   PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights);
433:   PetscFree2(Jinv, Jinvstar);
434:   return(0);
435: }

437: /*@C
438:   PetscQuadratureSetData - Sets the data defining the quadrature

440:   Not collective

442:   Input Parameters:
443: + q  - The PetscQuadrature object
444: . dim - The spatial dimension
445: . Nc - The number of components
446: . npoints - The number of quadrature points
447: . points - The coordinates of each quadrature point
448: - weights - The weight of each quadrature point

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

452:   Level: intermediate

454: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
455: @*/
456: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
457: {
460:   if (dim >= 0)     q->dim       = dim;
461:   if (Nc >= 0)      q->Nc        = Nc;
462:   if (npoints >= 0) q->numPoints = npoints;
463:   if (points) {
465:     q->points = points;
466:   }
467:   if (weights) {
469:     q->weights = weights;
470:   }
471:   return(0);
472: }

474: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
475: {
476:   PetscInt          q, d, c;
477:   PetscViewerFormat format;
478:   PetscErrorCode    ierr;

481:   if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D) with %D components\n", quad->order, quad->numPoints, quad->dim, quad->Nc);}
482:   else              {PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D)\n", quad->order, quad->numPoints, quad->dim);}
483:   PetscViewerGetFormat(v, &format);
484:   if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
485:   for (q = 0; q < quad->numPoints; ++q) {
486:     PetscViewerASCIIPrintf(v, "p%D (", q);
487:     PetscViewerASCIIUseTabs(v, PETSC_FALSE);
488:     for (d = 0; d < quad->dim; ++d) {
489:       if (d) {PetscViewerASCIIPrintf(v, ", ");}
490:       PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q*quad->dim+d]);
491:     }
492:     PetscViewerASCIIPrintf(v, ") ");
493:     if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "w%D (", q);}
494:     for (c = 0; c < quad->Nc; ++c) {
495:       if (c) {PetscViewerASCIIPrintf(v, ", ");}
496:       PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q*quad->Nc+c]);
497:     }
498:     if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, ")");}
499:     PetscViewerASCIIPrintf(v, "\n");
500:     PetscViewerASCIIUseTabs(v, PETSC_TRUE);
501:   }
502:   return(0);
503: }

505: /*@C
506:   PetscQuadratureView - Views a PetscQuadrature object

508:   Collective on quad

510:   Input Parameters:
511: + quad  - The PetscQuadrature object
512: - viewer - The PetscViewer object

514:   Level: beginner

516: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
517: @*/
518: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
519: {
520:   PetscBool      iascii;

526:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) quad), &viewer);}
527:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
528:   PetscViewerASCIIPushTab(viewer);
529:   if (iascii) {PetscQuadratureView_Ascii(quad, viewer);}
530:   PetscViewerASCIIPopTab(viewer);
531:   return(0);
532: }

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

537:   Not collective

539:   Input Parameter:
540: + q - The original PetscQuadrature
541: . numSubelements - The number of subelements the original element is divided into
542: . v0 - An array of the initial points for each subelement
543: - jac - An array of the Jacobian mappings from the reference to each subelement

545:   Output Parameters:
546: . dim - The dimension

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

550:  Not available from Fortran

552:   Level: intermediate

554: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
555: @*/
556: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
557: {
558:   const PetscReal *points,    *weights;
559:   PetscReal       *pointsRef, *weightsRef;
560:   PetscInt         dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
561:   PetscErrorCode   ierr;

568:   PetscQuadratureCreate(PETSC_COMM_SELF, qref);
569:   PetscQuadratureGetOrder(q, &order);
570:   PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
571:   npointsRef = npoints*numSubelements;
572:   PetscMalloc1(npointsRef*dim,&pointsRef);
573:   PetscMalloc1(npointsRef*Nc, &weightsRef);
574:   for (c = 0; c < numSubelements; ++c) {
575:     for (p = 0; p < npoints; ++p) {
576:       for (d = 0; d < dim; ++d) {
577:         pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
578:         for (e = 0; e < dim; ++e) {
579:           pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
580:         }
581:       }
582:       /* Could also use detJ here */
583:       for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
584:     }
585:   }
586:   PetscQuadratureSetOrder(*qref, order);
587:   PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
588:   return(0);
589: }

591: /*@
592:    PetscDTLegendreEval - evaluate Legendre polynomial at points

594:    Not Collective

596:    Input Arguments:
597: +  npoints - number of spatial points to evaluate at
598: .  points - array of locations to evaluate at
599: .  ndegree - number of basis degrees to evaluate
600: -  degrees - sorted array of degrees to evaluate

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

607:    Level: intermediate

609: .seealso: PetscDTGaussQuadrature()
610: @*/
611: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
612: {
613:   PetscInt i,maxdegree;

616:   if (!npoints || !ndegree) return(0);
617:   maxdegree = degrees[ndegree-1];
618:   for (i=0; i<npoints; i++) {
619:     PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
620:     PetscInt  j,k;
621:     x    = points[i];
622:     pm2  = 0;
623:     pm1  = 1;
624:     pd2  = 0;
625:     pd1  = 0;
626:     pdd2 = 0;
627:     pdd1 = 0;
628:     k    = 0;
629:     if (degrees[k] == 0) {
630:       if (B) B[i*ndegree+k] = pm1;
631:       if (D) D[i*ndegree+k] = pd1;
632:       if (D2) D2[i*ndegree+k] = pdd1;
633:       k++;
634:     }
635:     for (j=1; j<=maxdegree; j++,k++) {
636:       PetscReal p,d,dd;
637:       p    = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
638:       d    = pd2 + (2*j-1)*pm1;
639:       dd   = pdd2 + (2*j-1)*pd1;
640:       pm2  = pm1;
641:       pm1  = p;
642:       pd2  = pd1;
643:       pd1  = d;
644:       pdd2 = pdd1;
645:       pdd1 = dd;
646:       if (degrees[k] == j) {
647:         if (B) B[i*ndegree+k] = p;
648:         if (D) D[i*ndegree+k] = d;
649:         if (D2) D2[i*ndegree+k] = dd;
650:       }
651:     }
652:   }
653:   return(0);
654: }

656: /*@
657:    PetscDTGaussQuadrature - create Gauss quadrature

659:    Not Collective

661:    Input Arguments:
662: +  npoints - number of points
663: .  a - left end of interval (often-1)
664: -  b - right end of interval (often +1)

666:    Output Arguments:
667: +  x - quadrature points
668: -  w - quadrature weights

670:    Level: intermediate

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

675: .seealso: PetscDTLegendreEval()
676: @*/
677: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
678: {
680:   PetscInt       i;
681:   PetscReal      *work;
682:   PetscScalar    *Z;
683:   PetscBLASInt   N,LDZ,info;

686:   PetscCitationsRegister(GaussCitation, &GaussCite);
687:   /* Set up the Golub-Welsch system */
688:   for (i=0; i<npoints; i++) {
689:     x[i] = 0;                   /* diagonal is 0 */
690:     if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
691:   }
692:   PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
693:   PetscBLASIntCast(npoints,&N);
694:   LDZ  = N;
695:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
696:   PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
697:   PetscFPTrapPop();
698:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");

700:   for (i=0; i<(npoints+1)/2; i++) {
701:     PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
702:     x[i]           = (a+b)/2 - y*(b-a)/2;
703:     if (x[i] == -0.0) x[i] = 0.0;
704:     x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;

706:     w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
707:   }
708:   PetscFree2(Z,work);
709:   return(0);
710: }

712: static void qAndLEvaluation(PetscInt n, PetscReal x, PetscReal *q, PetscReal *qp, PetscReal *Ln)
713: /*
714:   Compute the polynomial q(x) = L_{N+1}(x) - L_{n-1}(x) and its derivative in
715:   addition to L_N(x) as these are needed for computing the GLL points via Newton's method.
716:   Reference: "Implementing Spectral Methods for Partial Differential Equations: Algorithms
717:   for Scientists and Engineers" by David A. Kopriva.
718: */
719: {
720:   PetscInt k;

722:   PetscReal Lnp;
723:   PetscReal Lnp1, Lnp1p;
724:   PetscReal Lnm1, Lnm1p;
725:   PetscReal Lnm2, Lnm2p;

727:   Lnm1  = 1.0;
728:   *Ln   = x;
729:   Lnm1p = 0.0;
730:   Lnp   = 1.0;

732:   for (k=2; k<=n; ++k) {
733:     Lnm2  = Lnm1;
734:     Lnm1  = *Ln;
735:     Lnm2p = Lnm1p;
736:     Lnm1p = Lnp;
737:     *Ln   = (2.*((PetscReal)k)-1.)/(1.0*((PetscReal)k))*x*Lnm1 - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm2;
738:     Lnp   = Lnm2p + (2.0*((PetscReal)k)-1.)*Lnm1;
739:   }
740:   k     = n+1;
741:   Lnp1  = (2.*((PetscReal)k)-1.)/(((PetscReal)k))*x*(*Ln) - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm1;
742:   Lnp1p = Lnm1p + (2.0*((PetscReal)k)-1.)*(*Ln);
743:   *q    = Lnp1 - Lnm1;
744:   *qp   = Lnp1p - Lnm1p;
745: }

747: /*@C
748:    PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre
749:                       nodes of a given size on the domain [-1,1]

751:    Not Collective

753:    Input Parameter:
754: +  n - number of grid nodes
755: -  type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON

757:    Output Arguments:
758: +  x - quadrature points
759: -  w - quadrature weights

761:    Notes:
762:     For n > 30  the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
763:           close enough to the desired solution

765:    These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes

767:    See  https://epubs.siam.org/doi/abs/10.1137/110855442  https://epubs.siam.org/doi/abs/10.1137/120889873 for better ways to compute GLL nodes

769:    Level: intermediate

771: .seealso: PetscDTGaussQuadrature()

773: @*/
774: PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal *x,PetscReal *w)
775: {

779:   if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide at least 2 grid points per element");

781:   if (type == PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA) {
782:     PetscReal      *M,si;
783:     PetscBLASInt   bn,lierr;
784:     PetscReal      x0,z0,z1,z2;
785:     PetscInt       i,p = npoints - 1,nn;

787:     x[0]   =-1.0;
788:     x[npoints-1] = 1.0;
789:     if (npoints-2 > 0){
790:       PetscMalloc1(npoints-1,&M);
791:       for (i=0; i<npoints-2; i++) {
792:         si  = ((PetscReal)i)+1.0;
793:         M[i]=0.5*PetscSqrtReal(si*(si+2.0)/((si+0.5)*(si+1.5)));
794:       }
795:       PetscBLASIntCast(npoints-2,&bn);
796:       PetscArrayzero(&x[1],bn);
797:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
798:       x0=0;
799:       PetscStackCallBLAS("LAPACKsteqr",LAPACKREALsteqr_("N",&bn,&x[1],M,&x0,&bn,M,&lierr));
800:       if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in STERF Lapack routine %d",(int)lierr);
801:       PetscFPTrapPop();
802:       PetscFree(M);
803:     }
804:     if ((npoints-1)%2==0) {
805:       x[(npoints-1)/2]   = 0.0; /* hard wire to exactly 0.0 since linear algebra produces nonzero */
806:     }

808:     w[0] = w[p] = 2.0/(((PetscReal)(p))*(((PetscReal)p)+1.0));
809:     z2 = -1.;                      /* Dummy value to avoid -Wmaybe-initialized */
810:     for (i=1; i<p; i++) {
811:       x0  = x[i];
812:       z0 = 1.0;
813:       z1 = x0;
814:       for (nn=1; nn<p; nn++) {
815:         z2 = x0*z1*(2.0*((PetscReal)nn)+1.0)/(((PetscReal)nn)+1.0)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.0));
816:         z0 = z1;
817:         z1 = z2;
818:       }
819:       w[i]=2.0/(((PetscReal)p)*(((PetscReal)p)+1.0)*z2*z2);
820:     }
821:   } else {
822:     PetscInt  j,m;
823:     PetscReal z1,z,q,qp,Ln;
824:     PetscReal *pt;
825:     PetscMalloc1(npoints,&pt);

827:     if (npoints > 30) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON produces incorrect answers for n > 30");
828:     x[0]     = -1.0;
829:     x[npoints-1]   = 1.0;
830:     w[0]   = w[npoints-1] = 2./(((PetscReal)npoints)*(((PetscReal)npoints)-1.0));
831:     m  = (npoints-1)/2; /* The roots are symmetric, so we only find half of them. */
832:     for (j=1; j<=m; j++) { /* Loop over the desired roots. */
833:       z = -1.0*PetscCosReal((PETSC_PI*((PetscReal)j)+0.25)/(((PetscReal)npoints)-1.0))-(3.0/(8.0*(((PetscReal)npoints)-1.0)*PETSC_PI))*(1.0/(((PetscReal)j)+0.25));
834:       /* Starting with the above approximation to the ith root, we enter */
835:       /* the main loop of refinement by Newton's method.                 */
836:       do {
837:         qAndLEvaluation(npoints-1,z,&q,&qp,&Ln);
838:         z1 = z;
839:         z  = z1-q/qp; /* Newton's method. */
840:       } while (PetscAbs(z-z1) > 10.*PETSC_MACHINE_EPSILON);
841:       qAndLEvaluation(npoints-1,z,&q,&qp,&Ln);

843:       x[j]       = z;
844:       x[npoints-1-j]   = -z;      /* and put in its symmetric counterpart.   */
845:       w[j]     = 2.0/(((PetscReal)npoints)*(((PetscReal)npoints)-1.)*Ln*Ln);  /* Compute the weight */
846:       w[npoints-1-j] = w[j];                 /* and its symmetric counterpart. */
847:       pt[j]=qp;
848:     }

850:     if ((npoints-1)%2==0) {
851:       qAndLEvaluation(npoints-1,0.0,&q,&qp,&Ln);
852:       x[(npoints-1)/2]   = 0.0;
853:       w[(npoints-1)/2] = 2.0/(((PetscReal)npoints)*(((PetscReal)npoints)-1.)*Ln*Ln);
854:     }
855:     PetscFree(pt);
856:   }
857:   return(0);
858: }

860: /*@
861:   PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature

863:   Not Collective

865:   Input Arguments:
866: + dim     - The spatial dimension
867: . Nc      - The number of components
868: . npoints - number of points in one dimension
869: . a       - left end of interval (often-1)
870: - b       - right end of interval (often +1)

872:   Output Argument:
873: . q - A PetscQuadrature object

875:   Level: intermediate

877: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
878: @*/
879: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
880: {
881:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
882:   PetscReal     *x, *w, *xw, *ww;

886:   PetscMalloc1(totpoints*dim,&x);
887:   PetscMalloc1(totpoints*Nc,&w);
888:   /* Set up the Golub-Welsch system */
889:   switch (dim) {
890:   case 0:
891:     PetscFree(x);
892:     PetscFree(w);
893:     PetscMalloc1(1, &x);
894:     PetscMalloc1(Nc, &w);
895:     x[0] = 0.0;
896:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
897:     break;
898:   case 1:
899:     PetscMalloc1(npoints,&ww);
900:     PetscDTGaussQuadrature(npoints, a, b, x, ww);
901:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
902:     PetscFree(ww);
903:     break;
904:   case 2:
905:     PetscMalloc2(npoints,&xw,npoints,&ww);
906:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
907:     for (i = 0; i < npoints; ++i) {
908:       for (j = 0; j < npoints; ++j) {
909:         x[(i*npoints+j)*dim+0] = xw[i];
910:         x[(i*npoints+j)*dim+1] = xw[j];
911:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
912:       }
913:     }
914:     PetscFree2(xw,ww);
915:     break;
916:   case 3:
917:     PetscMalloc2(npoints,&xw,npoints,&ww);
918:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
919:     for (i = 0; i < npoints; ++i) {
920:       for (j = 0; j < npoints; ++j) {
921:         for (k = 0; k < npoints; ++k) {
922:           x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
923:           x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
924:           x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
925:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
926:         }
927:       }
928:     }
929:     PetscFree2(xw,ww);
930:     break;
931:   default:
932:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
933:   }
934:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
935:   PetscQuadratureSetOrder(*q, 2*npoints-1);
936:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
937:   PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");
938:   return(0);
939: }

941: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
942:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
943: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
944: {
945:   PetscReal apb, pn1, pn2;
946:   PetscInt  k;

949:   if (!n) {*P = 1.0; return(0);}
950:   if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
951:   apb = a + b;
952:   pn2 = 1.0;
953:   pn1 = 0.5 * (a - b + (apb + 2.0) * x);
954:   *P  = 0.0;
955:   for (k = 2; k < n+1; ++k) {
956:     PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
957:     PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
958:     PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
959:     PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);

961:     a2  = a2 / a1;
962:     a3  = a3 / a1;
963:     a4  = a4 / a1;
964:     *P  = (a2 + a3 * x) * pn1 - a4 * pn2;
965:     pn2 = pn1;
966:     pn1 = *P;
967:   }
968:   return(0);
969: }

971: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
972: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
973: {
974:   PetscReal      nP;

978:   if (!n) {*P = 0.0; return(0);}
979:   PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
980:   *P   = 0.5 * (a + b + n + 1) * nP;
981:   return(0);
982: }

984: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
985: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
986: {
988:   *xi  = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
989:   *eta = y;
990:   return(0);
991: }

993: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
994: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
995: {
997:   *xi   = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
998:   *eta  = 0.5  * (1.0 + y) * (1.0 - z) - 1.0;
999:   *zeta = z;
1000:   return(0);
1001: }

1003: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
1004: {
1005:   PetscInt       maxIter = 100;
1006:   PetscReal      eps     = 1.0e-8;
1007:   PetscReal      a1, a2, a3, a4, a5, a6;
1008:   PetscInt       k;


1013:   a1      = PetscPowReal(2.0, a+b+1);
1014: #if defined(PETSC_HAVE_TGAMMA)
1015:   a2      = PetscTGamma(a + npoints + 1);
1016:   a3      = PetscTGamma(b + npoints + 1);
1017:   a4      = PetscTGamma(a + b + npoints + 1);
1018: #else
1019:   {
1020:     PetscInt ia, ib;

1022:     ia = (PetscInt) a;
1023:     ib = (PetscInt) b;
1024:     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 */
1025:       PetscDTFactorial(ia + npoints, &a2);
1026:       PetscDTFactorial(ib + npoints, &a3);
1027:       PetscDTFactorial(ia + ib + npoints, &a4);
1028:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
1029:   }
1030: #endif

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

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

1045:       for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
1046:       PetscDTComputeJacobi(a, b, npoints, r, &f);
1047:       PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
1048:       delta = f / (fp - f * s);
1049:       r     = r - delta;
1050:       if (PetscAbsReal(delta) < eps) break;
1051:     }
1052:     x[k] = r;
1053:     PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
1054:     w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
1055:   }
1056:   return(0);
1057: }

1059: /*@
1060:   PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex

1062:   Not Collective

1064:   Input Arguments:
1065: + dim     - The simplex dimension
1066: . Nc      - The number of components
1067: . npoints - The number of points in one dimension
1068: . a       - left end of interval (often-1)
1069: - b       - right end of interval (often +1)

1071:   Output Argument:
1072: . q - A PetscQuadrature object

1074:   Level: intermediate

1076:   References:
1077: .  1. - Karniadakis and Sherwin.  FIAT

1079: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
1080: @*/
1081: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1082: {
1083:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints;
1084:   PetscReal     *px, *wx, *py, *wy, *pz, *wz, *x, *w;
1085:   PetscInt       i, j, k, c;

1089:   if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
1090:   PetscMalloc1(totpoints*dim, &x);
1091:   PetscMalloc1(totpoints*Nc, &w);
1092:   switch (dim) {
1093:   case 0:
1094:     PetscFree(x);
1095:     PetscFree(w);
1096:     PetscMalloc1(1, &x);
1097:     PetscMalloc1(Nc, &w);
1098:     x[0] = 0.0;
1099:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
1100:     break;
1101:   case 1:
1102:     PetscMalloc1(npoints,&wx);
1103:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, wx);
1104:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i];
1105:     PetscFree(wx);
1106:     break;
1107:   case 2:
1108:     PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);
1109:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
1110:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
1111:     for (i = 0; i < npoints; ++i) {
1112:       for (j = 0; j < npoints; ++j) {
1113:         PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);
1114:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j];
1115:       }
1116:     }
1117:     PetscFree4(px,wx,py,wy);
1118:     break;
1119:   case 3:
1120:     PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);
1121:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
1122:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
1123:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);
1124:     for (i = 0; i < npoints; ++i) {
1125:       for (j = 0; j < npoints; ++j) {
1126:         for (k = 0; k < npoints; ++k) {
1127:           PetscDTMapCubeToTetrahedron_Internal(px[i], py[j], pz[k], &x[((i*npoints+j)*npoints+k)*3+0], &x[((i*npoints+j)*npoints+k)*3+1], &x[((i*npoints+j)*npoints+k)*3+2]);
1128:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k];
1129:         }
1130:       }
1131:     }
1132:     PetscFree6(px,wx,py,wy,pz,wz);
1133:     break;
1134:   default:
1135:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
1136:   }
1137:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
1138:   PetscQuadratureSetOrder(*q, 2*npoints-1);
1139:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
1140:   PetscObjectChangeTypeName((PetscObject)*q,"GaussJacobi");
1141:   return(0);
1142: }

1144: /*@
1145:   PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell

1147:   Not Collective

1149:   Input Arguments:
1150: + dim   - The cell dimension
1151: . level - The number of points in one dimension, 2^l
1152: . a     - left end of interval (often-1)
1153: - b     - right end of interval (often +1)

1155:   Output Argument:
1156: . q - A PetscQuadrature object

1158:   Level: intermediate

1160: .seealso: PetscDTGaussTensorQuadrature()
1161: @*/
1162: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
1163: {
1164:   const PetscInt  p     = 16;                        /* Digits of precision in the evaluation */
1165:   const PetscReal alpha = (b-a)/2.;                  /* Half-width of the integration interval */
1166:   const PetscReal beta  = (b+a)/2.;                  /* Center of the integration interval */
1167:   const PetscReal h     = PetscPowReal(2.0, -level); /* Step size, length between x_k */
1168:   PetscReal       xk;                                /* Quadrature point x_k on reference domain [-1, 1] */
1169:   PetscReal       wk    = 0.5*PETSC_PI;              /* Quadrature weight at x_k */
1170:   PetscReal      *x, *w;
1171:   PetscInt        K, k, npoints;
1172:   PetscErrorCode  ierr;

1175:   if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
1176:   if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
1177:   /* Find K such that the weights are < 32 digits of precision */
1178:   for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
1179:     wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
1180:   }
1181:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
1182:   PetscQuadratureSetOrder(*q, 2*K+1);
1183:   npoints = 2*K-1;
1184:   PetscMalloc1(npoints*dim, &x);
1185:   PetscMalloc1(npoints, &w);
1186:   /* Center term */
1187:   x[0] = beta;
1188:   w[0] = 0.5*alpha*PETSC_PI;
1189:   for (k = 1; k < K; ++k) {
1190:     wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1191:     xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
1192:     x[2*k-1] = -alpha*xk+beta;
1193:     w[2*k-1] = wk;
1194:     x[2*k+0] =  alpha*xk+beta;
1195:     w[2*k+0] = wk;
1196:   }
1197:   PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
1198:   return(0);
1199: }

1201: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1202: {
1203:   const PetscInt  p     = 16;        /* Digits of precision in the evaluation */
1204:   const PetscReal alpha = (b-a)/2.;  /* Half-width of the integration interval */
1205:   const PetscReal beta  = (b+a)/2.;  /* Center of the integration interval */
1206:   PetscReal       h     = 1.0;       /* Step size, length between x_k */
1207:   PetscInt        l     = 0;         /* Level of refinement, h = 2^{-l} */
1208:   PetscReal       osum  = 0.0;       /* Integral on last level */
1209:   PetscReal       psum  = 0.0;       /* Integral on the level before the last level */
1210:   PetscReal       sum;               /* Integral on current level */
1211:   PetscReal       yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1212:   PetscReal       lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1213:   PetscReal       wk;                /* Quadrature weight at x_k */
1214:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
1215:   PetscInt        d;                 /* Digits of precision in the integral */

1218:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1219:   /* Center term */
1220:   func(beta, &lval);
1221:   sum = 0.5*alpha*PETSC_PI*lval;
1222:   /* */
1223:   do {
1224:     PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
1225:     PetscInt  k = 1;

1227:     ++l;
1228:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1229:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1230:     psum = osum;
1231:     osum = sum;
1232:     h   *= 0.5;
1233:     sum *= 0.5;
1234:     do {
1235:       wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1236:       yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1237:       lx = -alpha*(1.0 - yk)+beta;
1238:       rx =  alpha*(1.0 - yk)+beta;
1239:       func(lx, &lval);
1240:       func(rx, &rval);
1241:       lterm   = alpha*wk*lval;
1242:       maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
1243:       sum    += lterm;
1244:       rterm   = alpha*wk*rval;
1245:       maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
1246:       sum    += rterm;
1247:       ++k;
1248:       /* Only need to evaluate every other point on refined levels */
1249:       if (l != 1) ++k;
1250:     } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */

1252:     d1 = PetscLog10Real(PetscAbsReal(sum - osum));
1253:     d2 = PetscLog10Real(PetscAbsReal(sum - psum));
1254:     d3 = PetscLog10Real(maxTerm) - p;
1255:     if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
1256:     else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
1257:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1258:   } while (d < digits && l < 12);
1259:   *sol = sum;

1261:   return(0);
1262: }

1264: #if defined(PETSC_HAVE_MPFR)
1265: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1266: {
1267:   const PetscInt  safetyFactor = 2;  /* Calculate abcissa until 2*p digits */
1268:   PetscInt        l            = 0;  /* Level of refinement, h = 2^{-l} */
1269:   mpfr_t          alpha;             /* Half-width of the integration interval */
1270:   mpfr_t          beta;              /* Center of the integration interval */
1271:   mpfr_t          h;                 /* Step size, length between x_k */
1272:   mpfr_t          osum;              /* Integral on last level */
1273:   mpfr_t          psum;              /* Integral on the level before the last level */
1274:   mpfr_t          sum;               /* Integral on current level */
1275:   mpfr_t          yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1276:   mpfr_t          lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1277:   mpfr_t          wk;                /* Quadrature weight at x_k */
1278:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
1279:   PetscInt        d;                 /* Digits of precision in the integral */
1280:   mpfr_t          pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;

1283:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1284:   /* Create high precision storage */
1285:   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);
1286:   /* Initialization */
1287:   mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
1288:   mpfr_set_d(beta,  0.5*(b+a), MPFR_RNDN);
1289:   mpfr_set_d(osum,  0.0,       MPFR_RNDN);
1290:   mpfr_set_d(psum,  0.0,       MPFR_RNDN);
1291:   mpfr_set_d(h,     1.0,       MPFR_RNDN);
1292:   mpfr_const_pi(pi2, MPFR_RNDN);
1293:   mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
1294:   /* Center term */
1295:   func(0.5*(b+a), &lval);
1296:   mpfr_set(sum, pi2, MPFR_RNDN);
1297:   mpfr_mul(sum, sum, alpha, MPFR_RNDN);
1298:   mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
1299:   /* */
1300:   do {
1301:     PetscReal d1, d2, d3, d4;
1302:     PetscInt  k = 1;

1304:     ++l;
1305:     mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
1306:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1307:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1308:     mpfr_set(psum, osum, MPFR_RNDN);
1309:     mpfr_set(osum,  sum, MPFR_RNDN);
1310:     mpfr_mul_d(h,   h,   0.5, MPFR_RNDN);
1311:     mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
1312:     do {
1313:       mpfr_set_si(kh, k, MPFR_RNDN);
1314:       mpfr_mul(kh, kh, h, MPFR_RNDN);
1315:       /* Weight */
1316:       mpfr_set(wk, h, MPFR_RNDN);
1317:       mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
1318:       mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
1319:       mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
1320:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
1321:       mpfr_sqr(tmp, tmp, MPFR_RNDN);
1322:       mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
1323:       mpfr_div(wk, wk, tmp, MPFR_RNDN);
1324:       /* Abscissa */
1325:       mpfr_set_d(yk, 1.0, MPFR_RNDZ);
1326:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
1327:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1328:       mpfr_exp(tmp, msinh, MPFR_RNDN);
1329:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1330:       /* Quadrature points */
1331:       mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
1332:       mpfr_mul(lx, lx, alpha, MPFR_RNDU);
1333:       mpfr_add(lx, lx, beta, MPFR_RNDU);
1334:       mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
1335:       mpfr_mul(rx, rx, alpha, MPFR_RNDD);
1336:       mpfr_add(rx, rx, beta, MPFR_RNDD);
1337:       /* Evaluation */
1338:       func(mpfr_get_d(lx, MPFR_RNDU), &lval);
1339:       func(mpfr_get_d(rx, MPFR_RNDD), &rval);
1340:       /* Update */
1341:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1342:       mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
1343:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
1344:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1345:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1346:       mpfr_set(curTerm, tmp, MPFR_RNDN);
1347:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1348:       mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
1349:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
1350:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1351:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1352:       mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
1353:       ++k;
1354:       /* Only need to evaluate every other point on refined levels */
1355:       if (l != 1) ++k;
1356:       mpfr_log10(tmp, wk, MPFR_RNDN);
1357:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1358:     } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
1359:     mpfr_sub(tmp, sum, osum, MPFR_RNDN);
1360:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1361:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1362:     d1 = mpfr_get_d(tmp, MPFR_RNDN);
1363:     mpfr_sub(tmp, sum, psum, MPFR_RNDN);
1364:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1365:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1366:     d2 = mpfr_get_d(tmp, MPFR_RNDN);
1367:     mpfr_log10(tmp, maxTerm, MPFR_RNDN);
1368:     d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
1369:     mpfr_log10(tmp, curTerm, MPFR_RNDN);
1370:     d4 = mpfr_get_d(tmp, MPFR_RNDN);
1371:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1372:   } while (d < digits && l < 8);
1373:   *sol = mpfr_get_d(sum, MPFR_RNDN);
1374:   /* Cleanup */
1375:   mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
1376:   return(0);
1377: }
1378: #else

1380: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1381: {
1382:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
1383: }
1384: #endif

1386: /* Overwrites A. Can only handle full-rank problems with m>=n
1387:  * A in column-major format
1388:  * Ainv in row-major format
1389:  * tau has length m
1390:  * worksize must be >= max(1,n)
1391:  */
1392: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1393: {
1395:   PetscBLASInt   M,N,K,lda,ldb,ldwork,info;
1396:   PetscScalar    *A,*Ainv,*R,*Q,Alpha;

1399: #if defined(PETSC_USE_COMPLEX)
1400:   {
1401:     PetscInt i,j;
1402:     PetscMalloc2(m*n,&A,m*n,&Ainv);
1403:     for (j=0; j<n; j++) {
1404:       for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1405:     }
1406:     mstride = m;
1407:   }
1408: #else
1409:   A = A_in;
1410:   Ainv = Ainv_out;
1411: #endif

1413:   PetscBLASIntCast(m,&M);
1414:   PetscBLASIntCast(n,&N);
1415:   PetscBLASIntCast(mstride,&lda);
1416:   PetscBLASIntCast(worksize,&ldwork);
1417:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1418:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1419:   PetscFPTrapPop();
1420:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1421:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1423:   /* Extract an explicit representation of Q */
1424:   Q = Ainv;
1425:   PetscArraycpy(Q,A,mstride*n);
1426:   K = N;                        /* full rank */
1427:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1428:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

1430:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1431:   Alpha = 1.0;
1432:   ldb = lda;
1433:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1434:   /* Ainv is Q, overwritten with inverse */

1436: #if defined(PETSC_USE_COMPLEX)
1437:   {
1438:     PetscInt i;
1439:     for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1440:     PetscFree2(A,Ainv);
1441:   }
1442: #endif
1443:   return(0);
1444: }

1446: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1447: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1448: {
1450:   PetscReal      *Bv;
1451:   PetscInt       i,j;

1454:   PetscMalloc1((ninterval+1)*ndegree,&Bv);
1455:   /* Point evaluation of L_p on all the source vertices */
1456:   PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1457:   /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1458:   for (i=0; i<ninterval; i++) {
1459:     for (j=0; j<ndegree; j++) {
1460:       if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1461:       else           B[i*ndegree+j]   = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1462:     }
1463:   }
1464:   PetscFree(Bv);
1465:   return(0);
1466: }

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

1471:    Not Collective

1473:    Input Arguments:
1474: +  degree - degree of reconstruction polynomial
1475: .  nsource - number of source intervals
1476: .  sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1477: .  ntarget - number of target intervals
1478: -  targetx - sorted coordinates of target cell boundaries (length ntarget+1)

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

1483:    Level: advanced

1485: .seealso: PetscDTLegendreEval()
1486: @*/
1487: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1488: {
1490:   PetscInt       i,j,k,*bdegrees,worksize;
1491:   PetscReal      xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1492:   PetscScalar    *tau,*work;

1498:   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);
1499: #if defined(PETSC_USE_DEBUG)
1500:   for (i=0; i<nsource; i++) {
1501:     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]);
1502:   }
1503:   for (i=0; i<ntarget; i++) {
1504:     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]);
1505:   }
1506: #endif
1507:   xmin = PetscMin(sourcex[0],targetx[0]);
1508:   xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1509:   center = (xmin + xmax)/2;
1510:   hscale = (xmax - xmin)/2;
1511:   worksize = nsource;
1512:   PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1513:   PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1514:   for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1515:   for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1516:   PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1517:   PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1518:   for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1519:   PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1520:   for (i=0; i<ntarget; i++) {
1521:     PetscReal rowsum = 0;
1522:     for (j=0; j<nsource; j++) {
1523:       PetscReal sum = 0;
1524:       for (k=0; k<degree+1; k++) {
1525:         sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1526:       }
1527:       R[i*nsource+j] = sum;
1528:       rowsum += sum;
1529:     }
1530:     for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1531:   }
1532:   PetscFree4(bdegrees,sourcey,Bsource,work);
1533:   PetscFree4(tau,Bsinv,targety,Btarget);
1534:   return(0);
1535: }

1537: /*@C
1538:    PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points

1540:    Not Collective

1542:    Input Parameter:
1543: +  n - the number of GLL nodes
1544: .  nodes - the GLL nodes
1545: .  weights - the GLL weights
1546: .  f - the function values at the nodes

1548:    Output Parameter:
1549: .  in - the value of the integral

1551:    Level: beginner

1553: .seealso: PetscDTGaussLobattoLegendreQuadrature()

1555: @*/
1556: PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal *nodes,PetscReal *weights,const PetscReal *f,PetscReal *in)
1557: {
1558:   PetscInt          i;

1561:   *in = 0.;
1562:   for (i=0; i<n; i++) {
1563:     *in += f[i]*f[i]*weights[i];
1564:   }
1565:   return(0);
1566: }

1568: /*@C
1569:    PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element

1571:    Not Collective

1573:    Input Parameter:
1574: +  n - the number of GLL nodes
1575: .  nodes - the GLL nodes
1576: .  weights - the GLL weights

1578:    Output Parameter:
1579: .  A - the stiffness element

1581:    Level: beginner

1583:    Notes:
1584:     Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy()

1586:    You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented (the array is symmetric)

1588: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()

1590: @*/
1591: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1592: {
1593:   PetscReal        **A;
1594:   PetscErrorCode  ierr;
1595:   const PetscReal  *gllnodes = nodes;
1596:   const PetscInt   p = n-1;
1597:   PetscReal        z0,z1,z2 = -1,x,Lpj,Lpr;
1598:   PetscInt         i,j,nn,r;

1601:   PetscMalloc1(n,&A);
1602:   PetscMalloc1(n*n,&A[0]);
1603:   for (i=1; i<n; i++) A[i] = A[i-1]+n;

1605:   for (j=1; j<p; j++) {
1606:     x  = gllnodes[j];
1607:     z0 = 1.;
1608:     z1 = x;
1609:     for (nn=1; nn<p; nn++) {
1610:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1611:       z0 = z1;
1612:       z1 = z2;
1613:     }
1614:     Lpj=z2;
1615:     for (r=1; r<p; r++) {
1616:       if (r == j) {
1617:         A[j][j]=2./(3.*(1.-gllnodes[j]*gllnodes[j])*Lpj*Lpj);
1618:       } else {
1619:         x  = gllnodes[r];
1620:         z0 = 1.;
1621:         z1 = x;
1622:         for (nn=1; nn<p; nn++) {
1623:           z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1624:           z0 = z1;
1625:           z1 = z2;
1626:         }
1627:         Lpr     = z2;
1628:         A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(gllnodes[j]-gllnodes[r])*(gllnodes[j]-gllnodes[r]));
1629:       }
1630:     }
1631:   }
1632:   for (j=1; j<p+1; j++) {
1633:     x  = gllnodes[j];
1634:     z0 = 1.;
1635:     z1 = x;
1636:     for (nn=1; nn<p; nn++) {
1637:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1638:       z0 = z1;
1639:       z1 = z2;
1640:     }
1641:     Lpj     = z2;
1642:     A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+gllnodes[j])*(1.+gllnodes[j]));
1643:     A[0][j] = A[j][0];
1644:   }
1645:   for (j=0; j<p; j++) {
1646:     x  = gllnodes[j];
1647:     z0 = 1.;
1648:     z1 = x;
1649:     for (nn=1; nn<p; nn++) {
1650:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1651:       z0 = z1;
1652:       z1 = z2;
1653:     }
1654:     Lpj=z2;

1656:     A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-gllnodes[j])*(1.-gllnodes[j]));
1657:     A[j][p] = A[p][j];
1658:   }
1659:   A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.;
1660:   A[p][p]=A[0][0];
1661:   *AA = A;
1662:   return(0);
1663: }

1665: /*@C
1666:    PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element

1668:    Not Collective

1670:    Input Parameter:
1671: +  n - the number of GLL nodes
1672: .  nodes - the GLL nodes
1673: .  weights - the GLL weightss
1674: -  A - the stiffness element

1676:    Level: beginner

1678: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate()

1680: @*/
1681: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1682: {

1686:   PetscFree((*AA)[0]);
1687:   PetscFree(*AA);
1688:   *AA  = NULL;
1689:   return(0);
1690: }

1692: /*@C
1693:    PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element

1695:    Not Collective

1697:    Input Parameter:
1698: +  n - the number of GLL nodes
1699: .  nodes - the GLL nodes
1700: .  weights - the GLL weights

1702:    Output Parameter:
1703: .  AA - the stiffness element
1704: -  AAT - the transpose of AA (pass in NULL if you do not need this array)

1706:    Level: beginner

1708:    Notes:
1709:     Destroy this with PetscGaussLobattoLegendreElementGradientDestroy()

1711:    You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented

1713: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()

1715: @*/
1716: PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
1717: {
1718:   PetscReal        **A, **AT = NULL;
1719:   PetscErrorCode  ierr;
1720:   const PetscReal  *gllnodes = nodes;
1721:   const PetscInt   p = n-1;
1722:   PetscReal        q,qp,Li, Lj,d0;
1723:   PetscInt         i,j;

1726:   PetscMalloc1(n,&A);
1727:   PetscMalloc1(n*n,&A[0]);
1728:   for (i=1; i<n; i++) A[i] = A[i-1]+n;

1730:   if (AAT) {
1731:     PetscMalloc1(n,&AT);
1732:     PetscMalloc1(n*n,&AT[0]);
1733:     for (i=1; i<n; i++) AT[i] = AT[i-1]+n;
1734:   }

1736:   if (n==1) {A[0][0] = 0.;}
1737:   d0 = (PetscReal)p*((PetscReal)p+1.)/4.;
1738:   for  (i=0; i<n; i++) {
1739:     for  (j=0; j<n; j++) {
1740:       A[i][j] = 0.;
1741:       qAndLEvaluation(p,gllnodes[i],&q,&qp,&Li);
1742:       qAndLEvaluation(p,gllnodes[j],&q,&qp,&Lj);
1743:       if (i!=j)             A[i][j] = Li/(Lj*(gllnodes[i]-gllnodes[j]));
1744:       if ((j==i) && (i==0)) A[i][j] = -d0;
1745:       if (j==i && i==p)     A[i][j] = d0;
1746:       if (AT) AT[j][i] = A[i][j];
1747:     }
1748:   }
1749:   if (AAT) *AAT = AT;
1750:   *AA  = A;
1751:   return(0);
1752: }

1754: /*@C
1755:    PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGaussLobattoLegendreElementGradientCreate()

1757:    Not Collective

1759:    Input Parameter:
1760: +  n - the number of GLL nodes
1761: .  nodes - the GLL nodes
1762: .  weights - the GLL weights
1763: .  AA - the stiffness element
1764: -  AAT - the transpose of the element

1766:    Level: beginner

1768: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionCreate()

1770: @*/
1771: PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
1772: {

1776:   PetscFree((*AA)[0]);
1777:   PetscFree(*AA);
1778:   *AA  = NULL;
1779:   if (*AAT) {
1780:     PetscFree((*AAT)[0]);
1781:     PetscFree(*AAT);
1782:     *AAT  = NULL;
1783:   }
1784:   return(0);
1785: }

1787: /*@C
1788:    PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element

1790:    Not Collective

1792:    Input Parameter:
1793: +  n - the number of GLL nodes
1794: .  nodes - the GLL nodes
1795: .  weights - the GLL weightss

1797:    Output Parameter:
1798: .  AA - the stiffness element

1800:    Level: beginner

1802:    Notes:
1803:     Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy()

1805:    This is the same as the Gradient operator multiplied by the diagonal mass matrix

1807:    You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented

1809: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionDestroy()

1811: @*/
1812: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1813: {
1814:   PetscReal       **D;
1815:   PetscErrorCode  ierr;
1816:   const PetscReal  *gllweights = weights;
1817:   const PetscInt   glln = n;
1818:   PetscInt         i,j;

1821:   PetscGaussLobattoLegendreElementGradientCreate(n,nodes,weights,&D,NULL);
1822:   for (i=0; i<glln; i++){
1823:     for (j=0; j<glln; j++) {
1824:       D[i][j] = gllweights[i]*D[i][j];
1825:     }
1826:   }
1827:   *AA = D;
1828:   return(0);
1829: }

1831: /*@C
1832:    PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element

1834:    Not Collective

1836:    Input Parameter:
1837: +  n - the number of GLL nodes
1838: .  nodes - the GLL nodes
1839: .  weights - the GLL weights
1840: -  A - advection

1842:    Level: beginner

1844: .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementAdvectionCreate()

1846: @*/
1847: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1848: {

1852:   PetscFree((*AA)[0]);
1853:   PetscFree(*AA);
1854:   *AA  = NULL;
1855:   return(0);
1856: }

1858: PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1859: {
1860:   PetscReal        **A;
1861:   PetscErrorCode  ierr;
1862:   const PetscReal  *gllweights = weights;
1863:   const PetscInt   glln = n;
1864:   PetscInt         i,j;

1867:   PetscMalloc1(glln,&A);
1868:   PetscMalloc1(glln*glln,&A[0]);
1869:   for (i=1; i<glln; i++) A[i] = A[i-1]+glln;
1870:   if (glln==1) {A[0][0] = 0.;}
1871:   for  (i=0; i<glln; i++) {
1872:     for  (j=0; j<glln; j++) {
1873:       A[i][j] = 0.;
1874:       if (j==i)     A[i][j] = gllweights[i];
1875:     }
1876:   }
1877:   *AA  = A;
1878:   return(0);
1879: }

1881: PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1882: {

1886:   PetscFree((*AA)[0]);
1887:   PetscFree(*AA);
1888:   *AA  = NULL;
1889:   return(0);
1890: }