Actual source code: dt.c

petsc-master 2019-08-19
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";

 25: /*@
 26:   PetscQuadratureCreate - Create a PetscQuadrature object

 28:   Collective

 30:   Input Parameter:
 31: . comm - The communicator for the PetscQuadrature object

 33:   Output Parameter:
 34: . q  - The PetscQuadrature object

 36:   Level: beginner

 38: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
 39: @*/
 40: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
 41: {

 46:   PetscSysInitializePackage();
 47:   PetscHeaderCreate(*q,PETSC_OBJECT_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
 48:   (*q)->dim       = -1;
 49:   (*q)->Nc        =  1;
 50:   (*q)->order     = -1;
 51:   (*q)->numPoints = 0;
 52:   (*q)->points    = NULL;
 53:   (*q)->weights   = NULL;
 54:   return(0);
 55: }

 57: /*@
 58:   PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object

 60:   Collective on q

 62:   Input Parameter:
 63: . q  - The PetscQuadrature object

 65:   Output Parameter:
 66: . r  - The new PetscQuadrature object

 68:   Level: beginner

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

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

 93: /*@
 94:   PetscQuadratureDestroy - Destroys a PetscQuadrature object

 96:   Collective on q

 98:   Input Parameter:
 99: . q  - The PetscQuadrature object

101:   Level: beginner

103: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
104: @*/
105: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
106: {

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

122: /*@
123:   PetscQuadratureGetOrder - Return the order of the method

125:   Not collective

127:   Input Parameter:
128: . q - The PetscQuadrature object

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

133:   Level: intermediate

135: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
136: @*/
137: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
138: {
142:   *order = q->order;
143:   return(0);
144: }

146: /*@
147:   PetscQuadratureSetOrder - Return the order of the method

149:   Not collective

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

155:   Level: intermediate

157: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
158: @*/
159: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
160: {
163:   q->order = order;
164:   return(0);
165: }

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

170:   Not collective

172:   Input Parameter:
173: . q - The PetscQuadrature object

175:   Output Parameter:
176: . Nc - The number of components

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

180:   Level: intermediate

182: .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
183: @*/
184: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
185: {
189:   *Nc = q->Nc;
190:   return(0);
191: }

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

196:   Not collective

198:   Input Parameters:
199: + q  - The PetscQuadrature object
200: - Nc - The number of components

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

204:   Level: intermediate

206: .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
207: @*/
208: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
209: {
212:   q->Nc = Nc;
213:   return(0);
214: }

216: /*@C
217:   PetscQuadratureGetData - Returns the data defining the quadrature

219:   Not collective

221:   Input Parameter:
222: . q  - The PetscQuadrature object

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

231:   Level: intermediate

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

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

265: /*@C
266:   PetscQuadratureSetData - Sets the data defining the quadrature

268:   Not collective

270:   Input Parameters:
271: + q  - The PetscQuadrature object
272: . dim - The spatial dimension
273: . Nc - The number of components
274: . npoints - The number of quadrature points
275: . points - The coordinates of each quadrature point
276: - weights - The weight of each quadrature point

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

280:   Level: intermediate

282: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
283: @*/
284: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
285: {
288:   if (dim >= 0)     q->dim       = dim;
289:   if (Nc >= 0)      q->Nc        = Nc;
290:   if (npoints >= 0) q->numPoints = npoints;
291:   if (points) {
293:     q->points = points;
294:   }
295:   if (weights) {
297:     q->weights = weights;
298:   }
299:   return(0);
300: }

302: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
303: {
304:   PetscInt          q, d, c;
305:   PetscViewerFormat format;
306:   PetscErrorCode    ierr;

309:   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);}
310:   else              {PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D)\n", quad->order, quad->numPoints, quad->dim);}
311:   PetscViewerGetFormat(v, &format);
312:   if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
313:   for (q = 0; q < quad->numPoints; ++q) {
314:     PetscViewerASCIIPrintf(v, "p%D (", q);
315:     PetscViewerASCIIUseTabs(v, PETSC_FALSE);
316:     for (d = 0; d < quad->dim; ++d) {
317:       if (d) {PetscViewerASCIIPrintf(v, ", ");}
318:       PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q*quad->dim+d]);
319:     }
320:     PetscViewerASCIIPrintf(v, ") ");
321:     if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "w%D (", q);}
322:     for (c = 0; c < quad->Nc; ++c) {
323:       if (c) {PetscViewerASCIIPrintf(v, ", ");}
324:       PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q*quad->Nc+c]);
325:     }
326:     if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, ")");}
327:     PetscViewerASCIIPrintf(v, "\n");
328:     PetscViewerASCIIUseTabs(v, PETSC_TRUE);
329:   }
330:   return(0);
331: }

333: /*@C
334:   PetscQuadratureView - Views a PetscQuadrature object

336:   Collective on quad

338:   Input Parameters:
339: + quad  - The PetscQuadrature object
340: - viewer - The PetscViewer object

342:   Level: beginner

344: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
345: @*/
346: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
347: {
348:   PetscBool      iascii;

354:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) quad), &viewer);}
355:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
356:   PetscViewerASCIIPushTab(viewer);
357:   if (iascii) {PetscQuadratureView_Ascii(quad, viewer);}
358:   PetscViewerASCIIPopTab(viewer);
359:   return(0);
360: }

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

365:   Not collective

367:   Input Parameter:
368: + q - The original PetscQuadrature
369: . numSubelements - The number of subelements the original element is divided into
370: . v0 - An array of the initial points for each subelement
371: - jac - An array of the Jacobian mappings from the reference to each subelement

373:   Output Parameters:
374: . dim - The dimension

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

378:  Not available from Fortran

380:   Level: intermediate

382: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
383: @*/
384: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
385: {
386:   const PetscReal *points,    *weights;
387:   PetscReal       *pointsRef, *weightsRef;
388:   PetscInt         dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
389:   PetscErrorCode   ierr;

396:   PetscQuadratureCreate(PETSC_COMM_SELF, qref);
397:   PetscQuadratureGetOrder(q, &order);
398:   PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
399:   npointsRef = npoints*numSubelements;
400:   PetscMalloc1(npointsRef*dim,&pointsRef);
401:   PetscMalloc1(npointsRef*Nc, &weightsRef);
402:   for (c = 0; c < numSubelements; ++c) {
403:     for (p = 0; p < npoints; ++p) {
404:       for (d = 0; d < dim; ++d) {
405:         pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
406:         for (e = 0; e < dim; ++e) {
407:           pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
408:         }
409:       }
410:       /* Could also use detJ here */
411:       for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
412:     }
413:   }
414:   PetscQuadratureSetOrder(*qref, order);
415:   PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
416:   return(0);
417: }

419: /*@
420:    PetscDTLegendreEval - evaluate Legendre polynomial at points

422:    Not Collective

424:    Input Arguments:
425: +  npoints - number of spatial points to evaluate at
426: .  points - array of locations to evaluate at
427: .  ndegree - number of basis degrees to evaluate
428: -  degrees - sorted array of degrees to evaluate

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

435:    Level: intermediate

437: .seealso: PetscDTGaussQuadrature()
438: @*/
439: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
440: {
441:   PetscInt i,maxdegree;

444:   if (!npoints || !ndegree) return(0);
445:   maxdegree = degrees[ndegree-1];
446:   for (i=0; i<npoints; i++) {
447:     PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
448:     PetscInt  j,k;
449:     x    = points[i];
450:     pm2  = 0;
451:     pm1  = 1;
452:     pd2  = 0;
453:     pd1  = 0;
454:     pdd2 = 0;
455:     pdd1 = 0;
456:     k    = 0;
457:     if (degrees[k] == 0) {
458:       if (B) B[i*ndegree+k] = pm1;
459:       if (D) D[i*ndegree+k] = pd1;
460:       if (D2) D2[i*ndegree+k] = pdd1;
461:       k++;
462:     }
463:     for (j=1; j<=maxdegree; j++,k++) {
464:       PetscReal p,d,dd;
465:       p    = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
466:       d    = pd2 + (2*j-1)*pm1;
467:       dd   = pdd2 + (2*j-1)*pd1;
468:       pm2  = pm1;
469:       pm1  = p;
470:       pd2  = pd1;
471:       pd1  = d;
472:       pdd2 = pdd1;
473:       pdd1 = dd;
474:       if (degrees[k] == j) {
475:         if (B) B[i*ndegree+k] = p;
476:         if (D) D[i*ndegree+k] = d;
477:         if (D2) D2[i*ndegree+k] = dd;
478:       }
479:     }
480:   }
481:   return(0);
482: }

484: /*@
485:    PetscDTGaussQuadrature - create Gauss quadrature

487:    Not Collective

489:    Input Arguments:
490: +  npoints - number of points
491: .  a - left end of interval (often-1)
492: -  b - right end of interval (often +1)

494:    Output Arguments:
495: +  x - quadrature points
496: -  w - quadrature weights

498:    Level: intermediate

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

503: .seealso: PetscDTLegendreEval()
504: @*/
505: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
506: {
508:   PetscInt       i;
509:   PetscReal      *work;
510:   PetscScalar    *Z;
511:   PetscBLASInt   N,LDZ,info;

514:   PetscCitationsRegister(GaussCitation, &GaussCite);
515:   /* Set up the Golub-Welsch system */
516:   for (i=0; i<npoints; i++) {
517:     x[i] = 0;                   /* diagonal is 0 */
518:     if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
519:   }
520:   PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
521:   PetscBLASIntCast(npoints,&N);
522:   LDZ  = N;
523:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
524:   PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
525:   PetscFPTrapPop();
526:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");

528:   for (i=0; i<(npoints+1)/2; i++) {
529:     PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
530:     x[i]           = (a+b)/2 - y*(b-a)/2;
531:     if (x[i] == -0.0) x[i] = 0.0;
532:     x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;

534:     w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
535:   }
536:   PetscFree2(Z,work);
537:   return(0);
538: }

540: static void qAndLEvaluation(PetscInt n, PetscReal x, PetscReal *q, PetscReal *qp, PetscReal *Ln)
541: /*
542:   Compute the polynomial q(x) = L_{N+1}(x) - L_{n-1}(x) and its derivative in
543:   addition to L_N(x) as these are needed for computing the GLL points via Newton's method.
544:   Reference: "Implementing Spectral Methods for Partial Differential Equations: Algorithms
545:   for Scientists and Engineers" by David A. Kopriva.
546: */
547: {
548:   PetscInt k;

550:   PetscReal Lnp;
551:   PetscReal Lnp1, Lnp1p;
552:   PetscReal Lnm1, Lnm1p;
553:   PetscReal Lnm2, Lnm2p;

555:   Lnm1  = 1.0;
556:   *Ln   = x;
557:   Lnm1p = 0.0;
558:   Lnp   = 1.0;

560:   for (k=2; k<=n; ++k) {
561:     Lnm2  = Lnm1;
562:     Lnm1  = *Ln;
563:     Lnm2p = Lnm1p;
564:     Lnm1p = Lnp;
565:     *Ln   = (2.*((PetscReal)k)-1.)/(1.0*((PetscReal)k))*x*Lnm1 - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm2;
566:     Lnp   = Lnm2p + (2.0*((PetscReal)k)-1.)*Lnm1;
567:   }
568:   k     = n+1;
569:   Lnp1  = (2.*((PetscReal)k)-1.)/(((PetscReal)k))*x*(*Ln) - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm1;
570:   Lnp1p = Lnm1p + (2.0*((PetscReal)k)-1.)*(*Ln);
571:   *q    = Lnp1 - Lnm1;
572:   *qp   = Lnp1p - Lnm1p;
573: }

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

579:    Not Collective

581:    Input Parameter:
582: +  n - number of grid nodes
583: -  type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON

585:    Output Arguments:
586: +  x - quadrature points
587: -  w - quadrature weights

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

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

595:    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

597:    Level: intermediate

599: .seealso: PetscDTGaussQuadrature()

601: @*/
602: PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal *x,PetscReal *w)
603: {

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

609:   if (type == PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA) {
610:     PetscReal      *M,si;
611:     PetscBLASInt   bn,lierr;
612:     PetscReal      x0,z0,z1,z2;
613:     PetscInt       i,p = npoints - 1,nn;

615:     x[0]   =-1.0;
616:     x[npoints-1] = 1.0;
617:     if (npoints-2 > 0){
618:       PetscMalloc1(npoints-1,&M);
619:       for (i=0; i<npoints-2; i++) {
620:         si  = ((PetscReal)i)+1.0;
621:         M[i]=0.5*PetscSqrtReal(si*(si+2.0)/((si+0.5)*(si+1.5)));
622:       }
623:       PetscBLASIntCast(npoints-2,&bn);
624:       PetscArrayzero(&x[1],bn);
625:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
626:       x0=0;
627:       PetscStackCallBLAS("LAPACKsteqr",LAPACKREALsteqr_("N",&bn,&x[1],M,&x0,&bn,M,&lierr));
628:       if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in STERF Lapack routine %d",(int)lierr);
629:       PetscFPTrapPop();
630:       PetscFree(M);
631:     }
632:     if ((npoints-1)%2==0) {
633:       x[(npoints-1)/2]   = 0.0; /* hard wire to exactly 0.0 since linear algebra produces nonzero */
634:     }

636:     w[0] = w[p] = 2.0/(((PetscReal)(p))*(((PetscReal)p)+1.0));
637:     z2 = -1.;                      /* Dummy value to avoid -Wmaybe-initialized */
638:     for (i=1; i<p; i++) {
639:       x0  = x[i];
640:       z0 = 1.0;
641:       z1 = x0;
642:       for (nn=1; nn<p; nn++) {
643:         z2 = x0*z1*(2.0*((PetscReal)nn)+1.0)/(((PetscReal)nn)+1.0)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.0));
644:         z0 = z1;
645:         z1 = z2;
646:       }
647:       w[i]=2.0/(((PetscReal)p)*(((PetscReal)p)+1.0)*z2*z2);
648:     }
649:   } else {
650:     PetscInt  j,m;
651:     PetscReal z1,z,q,qp,Ln;
652:     PetscReal *pt;
653:     PetscMalloc1(npoints,&pt);

655:     if (npoints > 30) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON produces incorrect answers for n > 30");
656:     x[0]     = -1.0;
657:     x[npoints-1]   = 1.0;
658:     w[0]   = w[npoints-1] = 2./(((PetscReal)npoints)*(((PetscReal)npoints)-1.0));;
659:     m  = (npoints-1)/2; /* The roots are symmetric, so we only find half of them. */
660:     for (j=1; j<=m; j++) { /* Loop over the desired roots. */
661:       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));
662:       /* Starting with the above approximation to the ith root, we enter */
663:       /* the main loop of refinement by Newton's method.                 */
664:       do {
665:         qAndLEvaluation(npoints-1,z,&q,&qp,&Ln);
666:         z1 = z;
667:         z  = z1-q/qp; /* Newton's method. */
668:       } while (PetscAbs(z-z1) > 10.*PETSC_MACHINE_EPSILON);
669:       qAndLEvaluation(npoints-1,z,&q,&qp,&Ln);

671:       x[j]       = z;
672:       x[npoints-1-j]   = -z;      /* and put in its symmetric counterpart.   */
673:       w[j]     = 2.0/(((PetscReal)npoints)*(((PetscReal)npoints)-1.)*Ln*Ln);  /* Compute the weight */
674:       w[npoints-1-j] = w[j];                 /* and its symmetric counterpart. */
675:       pt[j]=qp;
676:     }

678:     if ((npoints-1)%2==0) {
679:       qAndLEvaluation(npoints-1,0.0,&q,&qp,&Ln);
680:       x[(npoints-1)/2]   = 0.0;
681:       w[(npoints-1)/2] = 2.0/(((PetscReal)npoints)*(((PetscReal)npoints)-1.)*Ln*Ln);
682:     }
683:     PetscFree(pt);
684:   }
685:   return(0);
686: }

688: /*@
689:   PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature

691:   Not Collective

693:   Input Arguments:
694: + dim     - The spatial dimension
695: . Nc      - The number of components
696: . npoints - number of points in one dimension
697: . a       - left end of interval (often-1)
698: - b       - right end of interval (often +1)

700:   Output Argument:
701: . q - A PetscQuadrature object

703:   Level: intermediate

705: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
706: @*/
707: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
708: {
709:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
710:   PetscReal     *x, *w, *xw, *ww;

714:   PetscMalloc1(totpoints*dim,&x);
715:   PetscMalloc1(totpoints*Nc,&w);
716:   /* Set up the Golub-Welsch system */
717:   switch (dim) {
718:   case 0:
719:     PetscFree(x);
720:     PetscFree(w);
721:     PetscMalloc1(1, &x);
722:     PetscMalloc1(Nc, &w);
723:     x[0] = 0.0;
724:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
725:     break;
726:   case 1:
727:     PetscMalloc1(npoints,&ww);
728:     PetscDTGaussQuadrature(npoints, a, b, x, ww);
729:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
730:     PetscFree(ww);
731:     break;
732:   case 2:
733:     PetscMalloc2(npoints,&xw,npoints,&ww);
734:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
735:     for (i = 0; i < npoints; ++i) {
736:       for (j = 0; j < npoints; ++j) {
737:         x[(i*npoints+j)*dim+0] = xw[i];
738:         x[(i*npoints+j)*dim+1] = xw[j];
739:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
740:       }
741:     }
742:     PetscFree2(xw,ww);
743:     break;
744:   case 3:
745:     PetscMalloc2(npoints,&xw,npoints,&ww);
746:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
747:     for (i = 0; i < npoints; ++i) {
748:       for (j = 0; j < npoints; ++j) {
749:         for (k = 0; k < npoints; ++k) {
750:           x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
751:           x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
752:           x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
753:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
754:         }
755:       }
756:     }
757:     PetscFree2(xw,ww);
758:     break;
759:   default:
760:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
761:   }
762:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
763:   PetscQuadratureSetOrder(*q, 2*npoints-1);
764:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
765:   PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");
766:   return(0);
767: }

769: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
770:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
771: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
772: {
773:   PetscReal f = 1.0;
774:   PetscInt  i;

777:   for (i = 1; i < n+1; ++i) f *= i;
778:   *factorial = f;
779:   return(0);
780: }

782: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
783:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
784: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
785: {
786:   PetscReal apb, pn1, pn2;
787:   PetscInt  k;

790:   if (!n) {*P = 1.0; return(0);}
791:   if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
792:   apb = a + b;
793:   pn2 = 1.0;
794:   pn1 = 0.5 * (a - b + (apb + 2.0) * x);
795:   *P  = 0.0;
796:   for (k = 2; k < n+1; ++k) {
797:     PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
798:     PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
799:     PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
800:     PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);

802:     a2  = a2 / a1;
803:     a3  = a3 / a1;
804:     a4  = a4 / a1;
805:     *P  = (a2 + a3 * x) * pn1 - a4 * pn2;
806:     pn2 = pn1;
807:     pn1 = *P;
808:   }
809:   return(0);
810: }

812: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
813: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
814: {
815:   PetscReal      nP;

819:   if (!n) {*P = 0.0; return(0);}
820:   PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
821:   *P   = 0.5 * (a + b + n + 1) * nP;
822:   return(0);
823: }

825: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
826: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
827: {
829:   *xi  = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
830:   *eta = y;
831:   return(0);
832: }

834: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
835: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
836: {
838:   *xi   = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
839:   *eta  = 0.5  * (1.0 + y) * (1.0 - z) - 1.0;
840:   *zeta = z;
841:   return(0);
842: }

844: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
845: {
846:   PetscInt       maxIter = 100;
847:   PetscReal      eps     = 1.0e-8;
848:   PetscReal      a1, a2, a3, a4, a5, a6;
849:   PetscInt       k;


854:   a1      = PetscPowReal(2.0, a+b+1);
855: #if defined(PETSC_HAVE_TGAMMA)
856:   a2      = PetscTGamma(a + npoints + 1);
857:   a3      = PetscTGamma(b + npoints + 1);
858:   a4      = PetscTGamma(a + b + npoints + 1);
859: #else
860:   {
861:     PetscInt ia, ib;

863:     ia = (PetscInt) a;
864:     ib = (PetscInt) b;
865:     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 */
866:       PetscDTFactorial_Internal(ia + npoints, &a2);
867:       PetscDTFactorial_Internal(ib + npoints, &a3);
868:       PetscDTFactorial_Internal(ia + ib + npoints, &a4);
869:     } else {
870:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
871:     }
872:   }
873: #endif

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

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

888:       for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
889:       PetscDTComputeJacobi(a, b, npoints, r, &f);
890:       PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
891:       delta = f / (fp - f * s);
892:       r     = r - delta;
893:       if (PetscAbsReal(delta) < eps) break;
894:     }
895:     x[k] = r;
896:     PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
897:     w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
898:   }
899:   return(0);
900: }

902: /*@
903:   PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex

905:   Not Collective

907:   Input Arguments:
908: + dim     - The simplex dimension
909: . Nc      - The number of components
910: . npoints - The number of points in one dimension
911: . a       - left end of interval (often-1)
912: - b       - right end of interval (often +1)

914:   Output Argument:
915: . q - A PetscQuadrature object

917:   Level: intermediate

919:   References:
920: .  1. - Karniadakis and Sherwin.  FIAT

922: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
923: @*/
924: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
925: {
926:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints;
927:   PetscReal     *px, *wx, *py, *wy, *pz, *wz, *x, *w;
928:   PetscInt       i, j, k, c;

932:   if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
933:   PetscMalloc1(totpoints*dim, &x);
934:   PetscMalloc1(totpoints*Nc, &w);
935:   switch (dim) {
936:   case 0:
937:     PetscFree(x);
938:     PetscFree(w);
939:     PetscMalloc1(1, &x);
940:     PetscMalloc1(Nc, &w);
941:     x[0] = 0.0;
942:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
943:     break;
944:   case 1:
945:     PetscMalloc1(npoints,&wx);
946:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, wx);
947:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i];
948:     PetscFree(wx);
949:     break;
950:   case 2:
951:     PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);
952:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
953:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
954:     for (i = 0; i < npoints; ++i) {
955:       for (j = 0; j < npoints; ++j) {
956:         PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);
957:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j];
958:       }
959:     }
960:     PetscFree4(px,wx,py,wy);
961:     break;
962:   case 3:
963:     PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);
964:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
965:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
966:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);
967:     for (i = 0; i < npoints; ++i) {
968:       for (j = 0; j < npoints; ++j) {
969:         for (k = 0; k < npoints; ++k) {
970:           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]);
971:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k];
972:         }
973:       }
974:     }
975:     PetscFree6(px,wx,py,wy,pz,wz);
976:     break;
977:   default:
978:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
979:   }
980:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
981:   PetscQuadratureSetOrder(*q, 2*npoints-1);
982:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
983:   PetscObjectChangeTypeName((PetscObject)*q,"GaussJacobi");
984:   return(0);
985: }

987: /*@
988:   PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell

990:   Not Collective

992:   Input Arguments:
993: + dim   - The cell dimension
994: . level - The number of points in one dimension, 2^l
995: . a     - left end of interval (often-1)
996: - b     - right end of interval (often +1)

998:   Output Argument:
999: . q - A PetscQuadrature object

1001:   Level: intermediate

1003: .seealso: PetscDTGaussTensorQuadrature()
1004: @*/
1005: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
1006: {
1007:   const PetscInt  p     = 16;                        /* Digits of precision in the evaluation */
1008:   const PetscReal alpha = (b-a)/2.;                  /* Half-width of the integration interval */
1009:   const PetscReal beta  = (b+a)/2.;                  /* Center of the integration interval */
1010:   const PetscReal h     = PetscPowReal(2.0, -level); /* Step size, length between x_k */
1011:   PetscReal       xk;                                /* Quadrature point x_k on reference domain [-1, 1] */
1012:   PetscReal       wk    = 0.5*PETSC_PI;              /* Quadrature weight at x_k */
1013:   PetscReal      *x, *w;
1014:   PetscInt        K, k, npoints;
1015:   PetscErrorCode  ierr;

1018:   if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
1019:   if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
1020:   /* Find K such that the weights are < 32 digits of precision */
1021:   for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
1022:     wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
1023:   }
1024:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
1025:   PetscQuadratureSetOrder(*q, 2*K+1);
1026:   npoints = 2*K-1;
1027:   PetscMalloc1(npoints*dim, &x);
1028:   PetscMalloc1(npoints, &w);
1029:   /* Center term */
1030:   x[0] = beta;
1031:   w[0] = 0.5*alpha*PETSC_PI;
1032:   for (k = 1; k < K; ++k) {
1033:     wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1034:     xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
1035:     x[2*k-1] = -alpha*xk+beta;
1036:     w[2*k-1] = wk;
1037:     x[2*k+0] =  alpha*xk+beta;
1038:     w[2*k+0] = wk;
1039:   }
1040:   PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
1041:   return(0);
1042: }

1044: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1045: {
1046:   const PetscInt  p     = 16;        /* Digits of precision in the evaluation */
1047:   const PetscReal alpha = (b-a)/2.;  /* Half-width of the integration interval */
1048:   const PetscReal beta  = (b+a)/2.;  /* Center of the integration interval */
1049:   PetscReal       h     = 1.0;       /* Step size, length between x_k */
1050:   PetscInt        l     = 0;         /* Level of refinement, h = 2^{-l} */
1051:   PetscReal       osum  = 0.0;       /* Integral on last level */
1052:   PetscReal       psum  = 0.0;       /* Integral on the level before the last level */
1053:   PetscReal       sum;               /* Integral on current level */
1054:   PetscReal       yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1055:   PetscReal       lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1056:   PetscReal       wk;                /* Quadrature weight at x_k */
1057:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
1058:   PetscInt        d;                 /* Digits of precision in the integral */

1061:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1062:   /* Center term */
1063:   func(beta, &lval);
1064:   sum = 0.5*alpha*PETSC_PI*lval;
1065:   /* */
1066:   do {
1067:     PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
1068:     PetscInt  k = 1;

1070:     ++l;
1071:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1072:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1073:     psum = osum;
1074:     osum = sum;
1075:     h   *= 0.5;
1076:     sum *= 0.5;
1077:     do {
1078:       wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1079:       yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1080:       lx = -alpha*(1.0 - yk)+beta;
1081:       rx =  alpha*(1.0 - yk)+beta;
1082:       func(lx, &lval);
1083:       func(rx, &rval);
1084:       lterm   = alpha*wk*lval;
1085:       maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
1086:       sum    += lterm;
1087:       rterm   = alpha*wk*rval;
1088:       maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
1089:       sum    += rterm;
1090:       ++k;
1091:       /* Only need to evaluate every other point on refined levels */
1092:       if (l != 1) ++k;
1093:     } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */

1095:     d1 = PetscLog10Real(PetscAbsReal(sum - osum));
1096:     d2 = PetscLog10Real(PetscAbsReal(sum - psum));
1097:     d3 = PetscLog10Real(maxTerm) - p;
1098:     if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
1099:     else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
1100:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1101:   } while (d < digits && l < 12);
1102:   *sol = sum;

1104:   return(0);
1105: }

1107: #if defined(PETSC_HAVE_MPFR)
1108: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1109: {
1110:   const PetscInt  safetyFactor = 2;  /* Calculate abcissa until 2*p digits */
1111:   PetscInt        l            = 0;  /* Level of refinement, h = 2^{-l} */
1112:   mpfr_t          alpha;             /* Half-width of the integration interval */
1113:   mpfr_t          beta;              /* Center of the integration interval */
1114:   mpfr_t          h;                 /* Step size, length between x_k */
1115:   mpfr_t          osum;              /* Integral on last level */
1116:   mpfr_t          psum;              /* Integral on the level before the last level */
1117:   mpfr_t          sum;               /* Integral on current level */
1118:   mpfr_t          yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1119:   mpfr_t          lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1120:   mpfr_t          wk;                /* Quadrature weight at x_k */
1121:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
1122:   PetscInt        d;                 /* Digits of precision in the integral */
1123:   mpfr_t          pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;

1126:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1127:   /* Create high precision storage */
1128:   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);
1129:   /* Initialization */
1130:   mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
1131:   mpfr_set_d(beta,  0.5*(b+a), MPFR_RNDN);
1132:   mpfr_set_d(osum,  0.0,       MPFR_RNDN);
1133:   mpfr_set_d(psum,  0.0,       MPFR_RNDN);
1134:   mpfr_set_d(h,     1.0,       MPFR_RNDN);
1135:   mpfr_const_pi(pi2, MPFR_RNDN);
1136:   mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
1137:   /* Center term */
1138:   func(0.5*(b+a), &lval);
1139:   mpfr_set(sum, pi2, MPFR_RNDN);
1140:   mpfr_mul(sum, sum, alpha, MPFR_RNDN);
1141:   mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
1142:   /* */
1143:   do {
1144:     PetscReal d1, d2, d3, d4;
1145:     PetscInt  k = 1;

1147:     ++l;
1148:     mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
1149:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1150:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1151:     mpfr_set(psum, osum, MPFR_RNDN);
1152:     mpfr_set(osum,  sum, MPFR_RNDN);
1153:     mpfr_mul_d(h,   h,   0.5, MPFR_RNDN);
1154:     mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
1155:     do {
1156:       mpfr_set_si(kh, k, MPFR_RNDN);
1157:       mpfr_mul(kh, kh, h, MPFR_RNDN);
1158:       /* Weight */
1159:       mpfr_set(wk, h, MPFR_RNDN);
1160:       mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
1161:       mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
1162:       mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
1163:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
1164:       mpfr_sqr(tmp, tmp, MPFR_RNDN);
1165:       mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
1166:       mpfr_div(wk, wk, tmp, MPFR_RNDN);
1167:       /* Abscissa */
1168:       mpfr_set_d(yk, 1.0, MPFR_RNDZ);
1169:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
1170:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1171:       mpfr_exp(tmp, msinh, MPFR_RNDN);
1172:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1173:       /* Quadrature points */
1174:       mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
1175:       mpfr_mul(lx, lx, alpha, MPFR_RNDU);
1176:       mpfr_add(lx, lx, beta, MPFR_RNDU);
1177:       mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
1178:       mpfr_mul(rx, rx, alpha, MPFR_RNDD);
1179:       mpfr_add(rx, rx, beta, MPFR_RNDD);
1180:       /* Evaluation */
1181:       func(mpfr_get_d(lx, MPFR_RNDU), &lval);
1182:       func(mpfr_get_d(rx, MPFR_RNDD), &rval);
1183:       /* Update */
1184:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1185:       mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
1186:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
1187:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1188:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1189:       mpfr_set(curTerm, tmp, MPFR_RNDN);
1190:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1191:       mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
1192:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
1193:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1194:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1195:       mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
1196:       ++k;
1197:       /* Only need to evaluate every other point on refined levels */
1198:       if (l != 1) ++k;
1199:       mpfr_log10(tmp, wk, MPFR_RNDN);
1200:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1201:     } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
1202:     mpfr_sub(tmp, sum, osum, MPFR_RNDN);
1203:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1204:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1205:     d1 = mpfr_get_d(tmp, MPFR_RNDN);
1206:     mpfr_sub(tmp, sum, psum, MPFR_RNDN);
1207:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1208:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1209:     d2 = mpfr_get_d(tmp, MPFR_RNDN);
1210:     mpfr_log10(tmp, maxTerm, MPFR_RNDN);
1211:     d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
1212:     mpfr_log10(tmp, curTerm, MPFR_RNDN);
1213:     d4 = mpfr_get_d(tmp, MPFR_RNDN);
1214:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1215:   } while (d < digits && l < 8);
1216:   *sol = mpfr_get_d(sum, MPFR_RNDN);
1217:   /* Cleanup */
1218:   mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
1219:   return(0);
1220: }
1221: #else

1223: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1224: {
1225:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
1226: }
1227: #endif

1229: /* Overwrites A. Can only handle full-rank problems with m>=n
1230:  * A in column-major format
1231:  * Ainv in row-major format
1232:  * tau has length m
1233:  * worksize must be >= max(1,n)
1234:  */
1235: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1236: {
1238:   PetscBLASInt   M,N,K,lda,ldb,ldwork,info;
1239:   PetscScalar    *A,*Ainv,*R,*Q,Alpha;

1242: #if defined(PETSC_USE_COMPLEX)
1243:   {
1244:     PetscInt i,j;
1245:     PetscMalloc2(m*n,&A,m*n,&Ainv);
1246:     for (j=0; j<n; j++) {
1247:       for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1248:     }
1249:     mstride = m;
1250:   }
1251: #else
1252:   A = A_in;
1253:   Ainv = Ainv_out;
1254: #endif

1256:   PetscBLASIntCast(m,&M);
1257:   PetscBLASIntCast(n,&N);
1258:   PetscBLASIntCast(mstride,&lda);
1259:   PetscBLASIntCast(worksize,&ldwork);
1260:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1261:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1262:   PetscFPTrapPop();
1263:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1264:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1266:   /* Extract an explicit representation of Q */
1267:   Q = Ainv;
1268:   PetscArraycpy(Q,A,mstride*n);
1269:   K = N;                        /* full rank */
1270:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1271:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

1273:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1274:   Alpha = 1.0;
1275:   ldb = lda;
1276:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1277:   /* Ainv is Q, overwritten with inverse */

1279: #if defined(PETSC_USE_COMPLEX)
1280:   {
1281:     PetscInt i;
1282:     for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1283:     PetscFree2(A,Ainv);
1284:   }
1285: #endif
1286:   return(0);
1287: }

1289: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1290: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1291: {
1293:   PetscReal      *Bv;
1294:   PetscInt       i,j;

1297:   PetscMalloc1((ninterval+1)*ndegree,&Bv);
1298:   /* Point evaluation of L_p on all the source vertices */
1299:   PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1300:   /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1301:   for (i=0; i<ninterval; i++) {
1302:     for (j=0; j<ndegree; j++) {
1303:       if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1304:       else           B[i*ndegree+j]   = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1305:     }
1306:   }
1307:   PetscFree(Bv);
1308:   return(0);
1309: }

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

1314:    Not Collective

1316:    Input Arguments:
1317: +  degree - degree of reconstruction polynomial
1318: .  nsource - number of source intervals
1319: .  sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1320: .  ntarget - number of target intervals
1321: -  targetx - sorted coordinates of target cell boundaries (length ntarget+1)

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

1326:    Level: advanced

1328: .seealso: PetscDTLegendreEval()
1329: @*/
1330: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1331: {
1333:   PetscInt       i,j,k,*bdegrees,worksize;
1334:   PetscReal      xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1335:   PetscScalar    *tau,*work;

1341:   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);
1342: #if defined(PETSC_USE_DEBUG)
1343:   for (i=0; i<nsource; i++) {
1344:     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]);
1345:   }
1346:   for (i=0; i<ntarget; i++) {
1347:     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]);
1348:   }
1349: #endif
1350:   xmin = PetscMin(sourcex[0],targetx[0]);
1351:   xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1352:   center = (xmin + xmax)/2;
1353:   hscale = (xmax - xmin)/2;
1354:   worksize = nsource;
1355:   PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1356:   PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1357:   for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1358:   for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1359:   PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1360:   PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1361:   for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1362:   PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1363:   for (i=0; i<ntarget; i++) {
1364:     PetscReal rowsum = 0;
1365:     for (j=0; j<nsource; j++) {
1366:       PetscReal sum = 0;
1367:       for (k=0; k<degree+1; k++) {
1368:         sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1369:       }
1370:       R[i*nsource+j] = sum;
1371:       rowsum += sum;
1372:     }
1373:     for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1374:   }
1375:   PetscFree4(bdegrees,sourcey,Bsource,work);
1376:   PetscFree4(tau,Bsinv,targety,Btarget);
1377:   return(0);
1378: }

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

1383:    Not Collective

1385:    Input Parameter:
1386: +  n - the number of GLL nodes
1387: .  nodes - the GLL nodes
1388: .  weights - the GLL weights
1389: .  f - the function values at the nodes

1391:    Output Parameter:
1392: .  in - the value of the integral

1394:    Level: beginner

1396: .seealso: PetscDTGaussLobattoLegendreQuadrature()

1398: @*/
1399: PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal *nodes,PetscReal *weights,const PetscReal *f,PetscReal *in)
1400: {
1401:   PetscInt          i;

1404:   *in = 0.;
1405:   for (i=0; i<n; i++) {
1406:     *in += f[i]*f[i]*weights[i];
1407:   }
1408:   return(0);
1409: }

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

1414:    Not Collective

1416:    Input Parameter:
1417: +  n - the number of GLL nodes
1418: .  nodes - the GLL nodes
1419: .  weights - the GLL weights

1421:    Output Parameter:
1422: .  A - the stiffness element

1424:    Level: beginner

1426:    Notes:
1427:     Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy()

1429:    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)

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

1433: @*/
1434: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1435: {
1436:   PetscReal        **A;
1437:   PetscErrorCode  ierr;
1438:   const PetscReal  *gllnodes = nodes;
1439:   const PetscInt   p = n-1;
1440:   PetscReal        z0,z1,z2 = -1,x,Lpj,Lpr;
1441:   PetscInt         i,j,nn,r;

1444:   PetscMalloc1(n,&A);
1445:   PetscMalloc1(n*n,&A[0]);
1446:   for (i=1; i<n; i++) A[i] = A[i-1]+n;

1448:   for (j=1; j<p; j++) {
1449:     x  = gllnodes[j];
1450:     z0 = 1.;
1451:     z1 = x;
1452:     for (nn=1; nn<p; nn++) {
1453:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1454:       z0 = z1;
1455:       z1 = z2;
1456:     }
1457:     Lpj=z2;
1458:     for (r=1; r<p; r++) {
1459:       if (r == j) {
1460:         A[j][j]=2./(3.*(1.-gllnodes[j]*gllnodes[j])*Lpj*Lpj);
1461:       } else {
1462:         x  = gllnodes[r];
1463:         z0 = 1.;
1464:         z1 = x;
1465:         for (nn=1; nn<p; nn++) {
1466:           z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1467:           z0 = z1;
1468:           z1 = z2;
1469:         }
1470:         Lpr     = z2;
1471:         A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(gllnodes[j]-gllnodes[r])*(gllnodes[j]-gllnodes[r]));
1472:       }
1473:     }
1474:   }
1475:   for (j=1; j<p+1; j++) {
1476:     x  = gllnodes[j];
1477:     z0 = 1.;
1478:     z1 = x;
1479:     for (nn=1; nn<p; nn++) {
1480:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1481:       z0 = z1;
1482:       z1 = z2;
1483:     }
1484:     Lpj     = z2;
1485:     A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+gllnodes[j])*(1.+gllnodes[j]));
1486:     A[0][j] = A[j][0];
1487:   }
1488:   for (j=0; j<p; j++) {
1489:     x  = gllnodes[j];
1490:     z0 = 1.;
1491:     z1 = x;
1492:     for (nn=1; nn<p; nn++) {
1493:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
1494:       z0 = z1;
1495:       z1 = z2;
1496:     }
1497:     Lpj=z2;

1499:     A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-gllnodes[j])*(1.-gllnodes[j]));
1500:     A[j][p] = A[p][j];
1501:   }
1502:   A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.;
1503:   A[p][p]=A[0][0];
1504:   *AA = A;
1505:   return(0);
1506: }

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

1511:    Not Collective

1513:    Input Parameter:
1514: +  n - the number of GLL nodes
1515: .  nodes - the GLL nodes
1516: .  weights - the GLL weightss
1517: -  A - the stiffness element

1519:    Level: beginner

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

1523: @*/
1524: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1525: {

1529:   PetscFree((*AA)[0]);
1530:   PetscFree(*AA);
1531:   *AA  = NULL;
1532:   return(0);
1533: }

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

1538:    Not Collective

1540:    Input Parameter:
1541: +  n - the number of GLL nodes
1542: .  nodes - the GLL nodes
1543: .  weights - the GLL weights

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

1549:    Level: beginner

1551:    Notes:
1552:     Destroy this with PetscGaussLobattoLegendreElementGradientDestroy()

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

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

1558: @*/
1559: PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
1560: {
1561:   PetscReal        **A, **AT = NULL;
1562:   PetscErrorCode  ierr;
1563:   const PetscReal  *gllnodes = nodes;
1564:   const PetscInt   p = n-1;
1565:   PetscReal        q,qp,Li, Lj,d0;
1566:   PetscInt         i,j;

1569:   PetscMalloc1(n,&A);
1570:   PetscMalloc1(n*n,&A[0]);
1571:   for (i=1; i<n; i++) A[i] = A[i-1]+n;

1573:   if (AAT) {
1574:     PetscMalloc1(n,&AT);
1575:     PetscMalloc1(n*n,&AT[0]);
1576:     for (i=1; i<n; i++) AT[i] = AT[i-1]+n;
1577:   }

1579:   if (n==1) {A[0][0] = 0.;}
1580:   d0 = (PetscReal)p*((PetscReal)p+1.)/4.;
1581:   for  (i=0; i<n; i++) {
1582:     for  (j=0; j<n; j++) {
1583:       A[i][j] = 0.;
1584:       qAndLEvaluation(p,gllnodes[i],&q,&qp,&Li);
1585:       qAndLEvaluation(p,gllnodes[j],&q,&qp,&Lj);
1586:       if (i!=j)             A[i][j] = Li/(Lj*(gllnodes[i]-gllnodes[j]));
1587:       if ((j==i) && (i==0)) A[i][j] = -d0;
1588:       if (j==i && i==p)     A[i][j] = d0;
1589:       if (AT) AT[j][i] = A[i][j];
1590:     }
1591:   }
1592:   if (AAT) *AAT = AT;
1593:   *AA  = A;
1594:   return(0);
1595: }

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

1600:    Not Collective

1602:    Input Parameter:
1603: +  n - the number of GLL nodes
1604: .  nodes - the GLL nodes
1605: .  weights - the GLL weights
1606: .  AA - the stiffness element
1607: -  AAT - the transpose of the element

1609:    Level: beginner

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

1613: @*/
1614: PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
1615: {

1619:   PetscFree((*AA)[0]);
1620:   PetscFree(*AA);
1621:   *AA  = NULL;
1622:   if (*AAT) {
1623:     PetscFree((*AAT)[0]);
1624:     PetscFree(*AAT);
1625:     *AAT  = NULL;
1626:   }
1627:   return(0);
1628: }

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

1633:    Not Collective

1635:    Input Parameter:
1636: +  n - the number of GLL nodes
1637: .  nodes - the GLL nodes
1638: .  weights - the GLL weightss

1640:    Output Parameter:
1641: .  AA - the stiffness element

1643:    Level: beginner

1645:    Notes:
1646:     Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy()

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

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

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

1654: @*/
1655: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1656: {
1657:   PetscReal       **D;
1658:   PetscErrorCode  ierr;
1659:   const PetscReal  *gllweights = weights;
1660:   const PetscInt   glln = n;
1661:   PetscInt         i,j;

1664:   PetscGaussLobattoLegendreElementGradientCreate(n,nodes,weights,&D,NULL);
1665:   for (i=0; i<glln; i++){
1666:     for (j=0; j<glln; j++) {
1667:       D[i][j] = gllweights[i]*D[i][j];
1668:     }
1669:   }
1670:   *AA = D;
1671:   return(0);
1672: }

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

1677:    Not Collective

1679:    Input Parameter:
1680: +  n - the number of GLL nodes
1681: .  nodes - the GLL nodes
1682: .  weights - the GLL weights
1683: -  A - advection

1685:    Level: beginner

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

1689: @*/
1690: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1691: {

1695:   PetscFree((*AA)[0]);
1696:   PetscFree(*AA);
1697:   *AA  = NULL;
1698:   return(0);
1699: }

1701: PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1702: {
1703:   PetscReal        **A;
1704:   PetscErrorCode  ierr;
1705:   const PetscReal  *gllweights = weights;
1706:   const PetscInt   glln = n;
1707:   PetscInt         i,j;

1710:   PetscMalloc1(glln,&A);
1711:   PetscMalloc1(glln*glln,&A[0]);
1712:   for (i=1; i<glln; i++) A[i] = A[i-1]+glln;
1713:   if (glln==1) {A[0][0] = 0.;}
1714:   for  (i=0; i<glln; i++) {
1715:     for  (j=0; j<glln; j++) {
1716:       A[i][j] = 0.;
1717:       if (j==i)     A[i][j] = gllweights[i];
1718:     }
1719:   }
1720:   *AA  = A;
1721:   return(0);
1722: }

1724: PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
1725: {

1729:   PetscFree((*AA)[0]);
1730:   PetscFree(*AA);
1731:   *AA  = NULL;
1732:   return(0);
1733: }