Actual source code: dt.c
petsc-master 2020-01-21
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: }