Actual source code: dt.c
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 <petsc/private/petscfeimpl.h>
8: #include <petscviewer.h>
9: #include <petscdmplex.h>
10: #include <petscdmshell.h>
12: #if defined(PETSC_HAVE_MPFR)
13: #include <mpfr.h>
14: #endif
16: const char *const PetscDTNodeTypes_shifted[] = {"default", "gaussjacobi", "equispaced", "tanhsinh", "PETSCDTNODES_", NULL};
17: const char *const *const PetscDTNodeTypes = PetscDTNodeTypes_shifted + 1;
19: const char *const PetscDTSimplexQuadratureTypes_shifted[] = {"default", "conic", "minsym", "PETSCDTSIMPLEXQUAD_", NULL};
20: const char *const *const PetscDTSimplexQuadratureTypes = PetscDTSimplexQuadratureTypes_shifted + 1;
22: static PetscBool GolubWelschCite = PETSC_FALSE;
23: const char GolubWelschCitation[] = "@article{GolubWelsch1969,\n"
24: " author = {Golub and Welsch},\n"
25: " title = {Calculation of Quadrature Rules},\n"
26: " journal = {Math. Comp.},\n"
27: " volume = {23},\n"
28: " number = {106},\n"
29: " pages = {221--230},\n"
30: " year = {1969}\n}\n";
32: /* Numerical tests in src/dm/dt/tests/ex1.c show that when computing the nodes and weights of Gauss-Jacobi
33: quadrature rules:
35: - in double precision, Newton's method and Golub & Welsch both work for moderate degrees (< 100),
36: - in single precision, Newton's method starts producing incorrect roots around n = 15, but
37: the weights from Golub & Welsch become a problem before then: they produces errors
38: in computing the Jacobi-polynomial Gram matrix around n = 6.
40: So we default to Newton's method (required fewer dependencies) */
41: PetscBool PetscDTGaussQuadratureNewton_Internal = PETSC_TRUE;
43: PetscClassId PETSCQUADRATURE_CLASSID = 0;
45: /*@
46: PetscQuadratureCreate - Create a `PetscQuadrature` object
48: Collective
50: Input Parameter:
51: . comm - The communicator for the `PetscQuadrature` object
53: Output Parameter:
54: . q - The `PetscQuadrature` object
56: Level: beginner
58: .seealso: `PetscQuadrature`, `Petscquadraturedestroy()`, `PetscQuadratureGetData()`
59: @*/
60: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
61: {
62: PetscFunctionBegin;
64: PetscCall(DMInitializePackage());
65: PetscCall(PetscHeaderCreate(*q, PETSCQUADRATURE_CLASSID, "PetscQuadrature", "Quadrature", "DT", comm, PetscQuadratureDestroy, PetscQuadratureView));
66: (*q)->ct = DM_POLYTOPE_UNKNOWN;
67: (*q)->dim = -1;
68: (*q)->Nc = 1;
69: (*q)->order = -1;
70: (*q)->numPoints = 0;
71: (*q)->points = NULL;
72: (*q)->weights = NULL;
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /*@
77: PetscQuadratureDuplicate - Create a deep copy of the `PetscQuadrature` object
79: Collective
81: Input Parameter:
82: . q - The `PetscQuadrature` object
84: Output Parameter:
85: . r - The new `PetscQuadrature` object
87: Level: beginner
89: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureDestroy()`, `PetscQuadratureGetData()`
90: @*/
91: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
92: {
93: DMPolytopeType ct;
94: PetscInt order, dim, Nc, Nq;
95: const PetscReal *points, *weights;
96: PetscReal *p, *w;
98: PetscFunctionBegin;
100: PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)q), r));
101: PetscCall(PetscQuadratureGetCellType(q, &ct));
102: PetscCall(PetscQuadratureSetCellType(*r, ct));
103: PetscCall(PetscQuadratureGetOrder(q, &order));
104: PetscCall(PetscQuadratureSetOrder(*r, order));
105: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights));
106: PetscCall(PetscMalloc1(Nq * dim, &p));
107: PetscCall(PetscMalloc1(Nq * Nc, &w));
108: PetscCall(PetscArraycpy(p, points, Nq * dim));
109: PetscCall(PetscArraycpy(w, weights, Nc * Nq));
110: PetscCall(PetscQuadratureSetData(*r, dim, Nc, Nq, p, w));
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: /*@
115: PetscQuadratureDestroy - Destroys a `PetscQuadrature` object
117: Collective
119: Input Parameter:
120: . q - The `PetscQuadrature` object
122: Level: beginner
124: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
125: @*/
126: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
127: {
128: PetscFunctionBegin;
129: if (!*q) PetscFunctionReturn(PETSC_SUCCESS);
131: if (--((PetscObject)(*q))->refct > 0) {
132: *q = NULL;
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
135: PetscCall(PetscFree((*q)->points));
136: PetscCall(PetscFree((*q)->weights));
137: PetscCall(PetscHeaderDestroy(q));
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /*@
142: PetscQuadratureGetCellType - Return the cell type of the integration domain
144: Not Collective
146: Input Parameter:
147: . q - The `PetscQuadrature` object
149: Output Parameter:
150: . ct - The cell type of the integration domain
152: Level: intermediate
154: .seealso: `PetscQuadrature`, `PetscQuadratureSetCellType()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
155: @*/
156: PetscErrorCode PetscQuadratureGetCellType(PetscQuadrature q, DMPolytopeType *ct)
157: {
158: PetscFunctionBegin;
161: *ct = q->ct;
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: /*@
166: PetscQuadratureSetCellType - Set the cell type of the integration domain
168: Not Collective
170: Input Parameters:
171: + q - The `PetscQuadrature` object
172: - ct - The cell type of the integration domain
174: Level: intermediate
176: .seealso: `PetscQuadrature`, `PetscQuadratureGetCellType()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
177: @*/
178: PetscErrorCode PetscQuadratureSetCellType(PetscQuadrature q, DMPolytopeType ct)
179: {
180: PetscFunctionBegin;
182: q->ct = ct;
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*@
187: PetscQuadratureGetOrder - Return the order of the method in the `PetscQuadrature`
189: Not Collective
191: Input Parameter:
192: . q - The `PetscQuadrature` object
194: Output Parameter:
195: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
197: Level: intermediate
199: .seealso: `PetscQuadrature`, `PetscQuadratureSetOrder()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
200: @*/
201: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
202: {
203: PetscFunctionBegin;
206: *order = q->order;
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: /*@
211: PetscQuadratureSetOrder - Set the order of the method in the `PetscQuadrature`
213: Not Collective
215: Input Parameters:
216: + q - The `PetscQuadrature` object
217: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
219: Level: intermediate
221: .seealso: `PetscQuadrature`, `PetscQuadratureGetOrder()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
222: @*/
223: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
224: {
225: PetscFunctionBegin;
227: q->order = order;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated
234: Not Collective
236: Input Parameter:
237: . q - The `PetscQuadrature` object
239: Output Parameter:
240: . Nc - The number of components
242: Level: intermediate
244: Note:
245: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
247: .seealso: `PetscQuadrature`, `PetscQuadratureSetNumComponents()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
248: @*/
249: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
250: {
251: PetscFunctionBegin;
254: *Nc = q->Nc;
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /*@
259: PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated
261: Not Collective
263: Input Parameters:
264: + q - The `PetscQuadrature` object
265: - Nc - The number of components
267: Level: intermediate
269: Note:
270: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
272: .seealso: `PetscQuadrature`, `PetscQuadratureGetNumComponents()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
273: @*/
274: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
275: {
276: PetscFunctionBegin;
278: q->Nc = Nc;
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: /*@C
283: PetscQuadratureGetData - Returns the data defining the `PetscQuadrature`
285: Not Collective
287: Input Parameter:
288: . q - The `PetscQuadrature` object
290: Output Parameters:
291: + dim - The spatial dimension
292: . Nc - The number of components
293: . npoints - The number of quadrature points
294: . points - The coordinates of each quadrature point
295: - weights - The weight of each quadrature point
297: Level: intermediate
299: Fortran Note:
300: From Fortran you must call `PetscQuadratureRestoreData()` when you are done with the data
302: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureSetData()`
303: @*/
304: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
305: {
306: PetscFunctionBegin;
308: if (dim) {
310: *dim = q->dim;
311: }
312: if (Nc) {
314: *Nc = q->Nc;
315: }
316: if (npoints) {
318: *npoints = q->numPoints;
319: }
320: if (points) {
322: *points = q->points;
323: }
324: if (weights) {
326: *weights = q->weights;
327: }
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: /*@
332: PetscQuadratureEqual - determine whether two quadratures are equivalent
334: Input Parameters:
335: + A - A `PetscQuadrature` object
336: - B - Another `PetscQuadrature` object
338: Output Parameter:
339: . equal - `PETSC_TRUE` if the quadratures are the same
341: Level: intermediate
343: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`
344: @*/
345: PetscErrorCode PetscQuadratureEqual(PetscQuadrature A, PetscQuadrature B, PetscBool *equal)
346: {
347: PetscFunctionBegin;
351: *equal = PETSC_FALSE;
352: if (A->ct != B->ct || A->dim != B->dim || A->Nc != B->Nc || A->order != B->order || A->numPoints != B->numPoints) PetscFunctionReturn(PETSC_SUCCESS);
353: for (PetscInt i = 0; i < A->numPoints * A->dim; i++) {
354: if (PetscAbsReal(A->points[i] - B->points[i]) > PETSC_SMALL) PetscFunctionReturn(PETSC_SUCCESS);
355: }
356: if (!A->weights && !B->weights) {
357: *equal = PETSC_TRUE;
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
360: if (A->weights && B->weights) {
361: for (PetscInt i = 0; i < A->numPoints; i++) {
362: if (PetscAbsReal(A->weights[i] - B->weights[i]) > PETSC_SMALL) PetscFunctionReturn(PETSC_SUCCESS);
363: }
364: *equal = PETSC_TRUE;
365: }
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[])
370: {
371: PetscScalar *Js, *Jinvs;
372: PetscInt i, j, k;
373: PetscBLASInt bm, bn, info;
375: PetscFunctionBegin;
376: if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
377: PetscCall(PetscBLASIntCast(m, &bm));
378: PetscCall(PetscBLASIntCast(n, &bn));
379: #if defined(PETSC_USE_COMPLEX)
380: PetscCall(PetscMalloc2(m * n, &Js, m * n, &Jinvs));
381: for (i = 0; i < m * n; i++) Js[i] = J[i];
382: #else
383: Js = (PetscReal *)J;
384: Jinvs = Jinv;
385: #endif
386: if (m == n) {
387: PetscBLASInt *pivots;
388: PetscScalar *W;
390: PetscCall(PetscMalloc2(m, &pivots, m, &W));
392: PetscCall(PetscArraycpy(Jinvs, Js, m * m));
393: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info));
394: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscInt_FMT, (PetscInt)info);
395: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &info));
396: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscInt_FMT, (PetscInt)info);
397: PetscCall(PetscFree2(pivots, W));
398: } else if (m < n) {
399: PetscScalar *JJT;
400: PetscBLASInt *pivots;
401: PetscScalar *W;
403: PetscCall(PetscMalloc1(m * m, &JJT));
404: PetscCall(PetscMalloc2(m, &pivots, m, &W));
405: for (i = 0; i < m; i++) {
406: for (j = 0; j < m; j++) {
407: PetscScalar val = 0.;
409: for (k = 0; k < n; k++) val += Js[i * n + k] * Js[j * n + k];
410: JJT[i * m + j] = val;
411: }
412: }
414: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info));
415: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscInt_FMT, (PetscInt)info);
416: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info));
417: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscInt_FMT, (PetscInt)info);
418: for (i = 0; i < n; i++) {
419: for (j = 0; j < m; j++) {
420: PetscScalar val = 0.;
422: for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j];
423: Jinvs[i * m + j] = val;
424: }
425: }
426: PetscCall(PetscFree2(pivots, W));
427: PetscCall(PetscFree(JJT));
428: } else {
429: PetscScalar *JTJ;
430: PetscBLASInt *pivots;
431: PetscScalar *W;
433: PetscCall(PetscMalloc1(n * n, &JTJ));
434: PetscCall(PetscMalloc2(n, &pivots, n, &W));
435: for (i = 0; i < n; i++) {
436: for (j = 0; j < n; j++) {
437: PetscScalar val = 0.;
439: for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j];
440: JTJ[i * n + j] = val;
441: }
442: }
444: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bn, pivots, &info));
445: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscInt_FMT, (PetscInt)info);
446: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info));
447: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscInt_FMT, (PetscInt)info);
448: for (i = 0; i < n; i++) {
449: for (j = 0; j < m; j++) {
450: PetscScalar val = 0.;
452: for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k];
453: Jinvs[i * m + j] = val;
454: }
455: }
456: PetscCall(PetscFree2(pivots, W));
457: PetscCall(PetscFree(JTJ));
458: }
459: #if defined(PETSC_USE_COMPLEX)
460: for (i = 0; i < m * n; i++) Jinv[i] = PetscRealPart(Jinvs[i]);
461: PetscCall(PetscFree2(Js, Jinvs));
462: #endif
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
466: /*@
467: PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation.
469: Collective
471: Input Parameters:
472: + q - the quadrature functional
473: . imageDim - the dimension of the image of the transformation
474: . origin - a point in the original space
475: . originImage - the image of the origin under the transformation
476: . J - the Jacobian of the image: an [imageDim x dim] matrix in row major order
477: - 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]
479: Output Parameter:
480: . 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.
482: Level: intermediate
484: Note:
485: 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.
487: .seealso: `PetscQuadrature`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
488: @*/
489: PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq)
490: {
491: PetscInt dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c;
492: const PetscReal *points;
493: const PetscReal *weights;
494: PetscReal *imagePoints, *imageWeights;
495: PetscReal *Jinv;
496: PetscReal *Jinvstar;
498: PetscFunctionBegin;
500: PetscCheck(imageDim >= PetscAbsInt(formDegree), PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Cannot represent a %" PetscInt_FMT "-form in %" PetscInt_FMT " dimensions", PetscAbsInt(formDegree), imageDim);
501: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights));
502: PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize));
503: PetscCheck(Nc % formSize == 0, PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Number of components %" PetscInt_FMT " is not a multiple of formSize %" PetscInt_FMT, Nc, formSize);
504: Ncopies = Nc / formSize;
505: PetscCall(PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize));
506: imageNc = Ncopies * imageFormSize;
507: PetscCall(PetscMalloc1(Npoints * imageDim, &imagePoints));
508: PetscCall(PetscMalloc1(Npoints * imageNc, &imageWeights));
509: PetscCall(PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar));
510: PetscCall(PetscDTJacobianInverse_Internal(imageDim, dim, J, Jinv));
511: PetscCall(PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar));
512: for (pt = 0; pt < Npoints; pt++) {
513: const PetscReal *point = &points[pt * dim];
514: PetscReal *imagePoint = &imagePoints[pt * imageDim];
516: for (i = 0; i < imageDim; i++) {
517: PetscReal val = originImage[i];
519: for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]);
520: imagePoint[i] = val;
521: }
522: for (c = 0; c < Ncopies; c++) {
523: const PetscReal *form = &weights[pt * Nc + c * formSize];
524: PetscReal *imageForm = &imageWeights[pt * imageNc + c * imageFormSize];
526: for (i = 0; i < imageFormSize; i++) {
527: PetscReal val = 0.;
529: for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j];
530: imageForm[i] = val;
531: }
532: }
533: }
534: PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq));
535: PetscCall(PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights));
536: PetscCall(PetscFree2(Jinv, Jinvstar));
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: /*@C
541: PetscQuadratureSetData - Sets the data defining the quadrature
543: Not Collective
545: Input Parameters:
546: + q - The `PetscQuadrature` object
547: . dim - The spatial dimension
548: . Nc - The number of components
549: . npoints - The number of quadrature points
550: . points - The coordinates of each quadrature point
551: - weights - The weight of each quadrature point
553: Level: intermediate
555: Note:
556: This routine owns the references to points and weights, so they must be allocated using `PetscMalloc()` and the user should not free them.
558: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
559: @*/
560: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
561: {
562: PetscFunctionBegin;
564: if (dim >= 0) q->dim = dim;
565: if (Nc >= 0) q->Nc = Nc;
566: if (npoints >= 0) q->numPoints = npoints;
567: if (points) {
569: q->points = points;
570: }
571: if (weights) {
573: q->weights = weights;
574: }
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
579: {
580: PetscInt q, d, c;
581: PetscViewerFormat format;
583: PetscFunctionBegin;
584: if (quad->Nc > 1)
585: PetscCall(PetscViewerASCIIPrintf(v, "Quadrature on a %s of order %" PetscInt_FMT " on %" PetscInt_FMT " points (dim %" PetscInt_FMT ") with %" PetscInt_FMT " components\n", DMPolytopeTypes[quad->ct], quad->order, quad->numPoints, quad->dim, quad->Nc));
586: else PetscCall(PetscViewerASCIIPrintf(v, "Quadrature on a %s of order %" PetscInt_FMT " on %" PetscInt_FMT " points (dim %" PetscInt_FMT ")\n", DMPolytopeTypes[quad->ct], quad->order, quad->numPoints, quad->dim));
587: PetscCall(PetscViewerGetFormat(v, &format));
588: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
589: for (q = 0; q < quad->numPoints; ++q) {
590: PetscCall(PetscViewerASCIIPrintf(v, "p%" PetscInt_FMT " (", q));
591: PetscCall(PetscViewerASCIIUseTabs(v, PETSC_FALSE));
592: for (d = 0; d < quad->dim; ++d) {
593: if (d) PetscCall(PetscViewerASCIIPrintf(v, ", "));
594: PetscCall(PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q * quad->dim + d]));
595: }
596: PetscCall(PetscViewerASCIIPrintf(v, ") "));
597: if (quad->Nc > 1) PetscCall(PetscViewerASCIIPrintf(v, "w%" PetscInt_FMT " (", q));
598: for (c = 0; c < quad->Nc; ++c) {
599: if (c) PetscCall(PetscViewerASCIIPrintf(v, ", "));
600: PetscCall(PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q * quad->Nc + c]));
601: }
602: if (quad->Nc > 1) PetscCall(PetscViewerASCIIPrintf(v, ")"));
603: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
604: PetscCall(PetscViewerASCIIUseTabs(v, PETSC_TRUE));
605: }
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: /*@C
610: PetscQuadratureView - View a `PetscQuadrature` object
612: Collective
614: Input Parameters:
615: + quad - The `PetscQuadrature` object
616: - viewer - The `PetscViewer` object
618: Level: beginner
620: .seealso: `PetscQuadrature`, `PetscViewer`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
621: @*/
622: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
623: {
624: PetscBool iascii;
626: PetscFunctionBegin;
629: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)quad), &viewer));
630: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
631: PetscCall(PetscViewerASCIIPushTab(viewer));
632: if (iascii) PetscCall(PetscQuadratureView_Ascii(quad, viewer));
633: PetscCall(PetscViewerASCIIPopTab(viewer));
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
637: /*@C
638: PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement
640: Not Collective; No Fortran Support
642: Input Parameters:
643: + q - The original `PetscQuadrature`
644: . numSubelements - The number of subelements the original element is divided into
645: . v0 - An array of the initial points for each subelement
646: - jac - An array of the Jacobian mappings from the reference to each subelement
648: Output Parameter:
649: . dim - The dimension
651: Level: intermediate
653: Note:
654: Together v0 and jac define an affine mapping from the original reference element to each subelement
656: .seealso: `PetscQuadrature`, `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()`
657: @*/
658: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
659: {
660: DMPolytopeType ct;
661: const PetscReal *points, *weights;
662: PetscReal *pointsRef, *weightsRef;
663: PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
665: PetscFunctionBegin;
670: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, qref));
671: PetscCall(PetscQuadratureGetCellType(q, &ct));
672: PetscCall(PetscQuadratureGetOrder(q, &order));
673: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
674: npointsRef = npoints * numSubelements;
675: PetscCall(PetscMalloc1(npointsRef * dim, &pointsRef));
676: PetscCall(PetscMalloc1(npointsRef * Nc, &weightsRef));
677: for (c = 0; c < numSubelements; ++c) {
678: for (p = 0; p < npoints; ++p) {
679: for (d = 0; d < dim; ++d) {
680: pointsRef[(c * npoints + p) * dim + d] = v0[c * dim + d];
681: for (e = 0; e < dim; ++e) pointsRef[(c * npoints + p) * dim + d] += jac[(c * dim + d) * dim + e] * (points[p * dim + e] + 1.0);
682: }
683: /* Could also use detJ here */
684: for (cp = 0; cp < Nc; ++cp) weightsRef[(c * npoints + p) * Nc + cp] = weights[p * Nc + cp] / numSubelements;
685: }
686: }
687: PetscCall(PetscQuadratureSetCellType(*qref, ct));
688: PetscCall(PetscQuadratureSetOrder(*qref, order));
689: PetscCall(PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef));
690: PetscFunctionReturn(PETSC_SUCCESS);
691: }
693: /* Compute the coefficients for the Jacobi polynomial recurrence,
694: *
695: * J^{a,b}_n(x) = (cnm1 + cnm1x * x) * J^{a,b}_{n-1}(x) - cnm2 * J^{a,b}_{n-2}(x).
696: */
697: #define PetscDTJacobiRecurrence_Internal(n, a, b, cnm1, cnm1x, cnm2) \
698: do { \
699: PetscReal _a = (a); \
700: PetscReal _b = (b); \
701: PetscReal _n = (n); \
702: if (n == 1) { \
703: (cnm1) = (_a - _b) * 0.5; \
704: (cnm1x) = (_a + _b + 2.) * 0.5; \
705: (cnm2) = 0.; \
706: } else { \
707: PetscReal _2n = _n + _n; \
708: PetscReal _d = (_2n * (_n + _a + _b) * (_2n + _a + _b - 2)); \
709: PetscReal _n1 = (_2n + _a + _b - 1.) * (_a * _a - _b * _b); \
710: PetscReal _n1x = (_2n + _a + _b - 1.) * (_2n + _a + _b) * (_2n + _a + _b - 2); \
711: PetscReal _n2 = 2. * ((_n + _a - 1.) * (_n + _b - 1.) * (_2n + _a + _b)); \
712: (cnm1) = _n1 / _d; \
713: (cnm1x) = _n1x / _d; \
714: (cnm2) = _n2 / _d; \
715: } \
716: } while (0)
718: /*@
719: PetscDTJacobiNorm - Compute the weighted L2 norm of a Jacobi polynomial.
721: $\| P^{\alpha,\beta}_n \|_{\alpha,\beta}^2 = \int_{-1}^1 (1 + x)^{\alpha} (1 - x)^{\beta} P^{\alpha,\beta}_n (x)^2 dx.$
723: Input Parameters:
724: - alpha - the left exponent > -1
725: . beta - the right exponent > -1
726: + n - the polynomial degree
728: Output Parameter:
729: . norm - the weighted L2 norm
731: Level: beginner
733: .seealso: `PetscQuadrature`, `PetscDTJacobiEval()`
734: @*/
735: PetscErrorCode PetscDTJacobiNorm(PetscReal alpha, PetscReal beta, PetscInt n, PetscReal *norm)
736: {
737: PetscReal twoab1;
738: PetscReal gr;
740: PetscFunctionBegin;
741: PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent alpha %g <= -1. invalid", (double)alpha);
742: PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent beta %g <= -1. invalid", (double)beta);
743: PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "n %" PetscInt_FMT " < 0 invalid", n);
744: twoab1 = PetscPowReal(2., alpha + beta + 1.);
745: #if defined(PETSC_HAVE_LGAMMA)
746: if (!n) {
747: gr = PetscExpReal(PetscLGamma(alpha + 1.) + PetscLGamma(beta + 1.) - PetscLGamma(alpha + beta + 2.));
748: } else {
749: gr = PetscExpReal(PetscLGamma(n + alpha + 1.) + PetscLGamma(n + beta + 1.) - (PetscLGamma(n + 1.) + PetscLGamma(n + alpha + beta + 1.))) / (n + n + alpha + beta + 1.);
750: }
751: #else
752: {
753: PetscInt alphai = (PetscInt)alpha;
754: PetscInt betai = (PetscInt)beta;
755: PetscInt i;
757: gr = n ? (1. / (n + n + alpha + beta + 1.)) : 1.;
758: if ((PetscReal)alphai == alpha) {
759: if (!n) {
760: for (i = 0; i < alphai; i++) gr *= (i + 1.) / (beta + i + 1.);
761: gr /= (alpha + beta + 1.);
762: } else {
763: for (i = 0; i < alphai; i++) gr *= (n + i + 1.) / (n + beta + i + 1.);
764: }
765: } else if ((PetscReal)betai == beta) {
766: if (!n) {
767: for (i = 0; i < betai; i++) gr *= (i + 1.) / (alpha + i + 2.);
768: gr /= (alpha + beta + 1.);
769: } else {
770: for (i = 0; i < betai; i++) gr *= (n + i + 1.) / (n + alpha + i + 1.);
771: }
772: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
773: }
774: #endif
775: *norm = PetscSqrtReal(twoab1 * gr);
776: PetscFunctionReturn(PETSC_SUCCESS);
777: }
779: static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p)
780: {
781: PetscReal ak, bk;
782: PetscReal abk1;
783: PetscInt i, l, maxdegree;
785: PetscFunctionBegin;
786: maxdegree = degrees[ndegree - 1] - k;
787: ak = a + k;
788: bk = b + k;
789: abk1 = a + b + k + 1.;
790: if (maxdegree < 0) {
791: for (i = 0; i < npoints; i++)
792: for (l = 0; l < ndegree; l++) p[i * ndegree + l] = 0.;
793: PetscFunctionReturn(PETSC_SUCCESS);
794: }
795: for (i = 0; i < npoints; i++) {
796: PetscReal pm1, pm2, x;
797: PetscReal cnm1, cnm1x, cnm2;
798: PetscInt j, m;
800: x = points[i];
801: pm2 = 1.;
802: PetscDTJacobiRecurrence_Internal(1, ak, bk, cnm1, cnm1x, cnm2);
803: pm1 = (cnm1 + cnm1x * x);
804: l = 0;
805: while (l < ndegree && degrees[l] - k < 0) p[l++] = 0.;
806: while (l < ndegree && degrees[l] - k == 0) {
807: p[l] = pm2;
808: for (m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5;
809: l++;
810: }
811: while (l < ndegree && degrees[l] - k == 1) {
812: p[l] = pm1;
813: for (m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5;
814: l++;
815: }
816: for (j = 2; j <= maxdegree; j++) {
817: PetscReal pp;
819: PetscDTJacobiRecurrence_Internal(j, ak, bk, cnm1, cnm1x, cnm2);
820: pp = (cnm1 + cnm1x * x) * pm1 - cnm2 * pm2;
821: pm2 = pm1;
822: pm1 = pp;
823: while (l < ndegree && degrees[l] - k == j) {
824: p[l] = pp;
825: for (m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5;
826: l++;
827: }
828: }
829: p += ndegree;
830: }
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: /*@
835: PetscDTJacobiEvalJet - Evaluate the jet (function and derivatives) of the Jacobi polynomials basis up to a given degree.
836: The Jacobi polynomials with indices $\alpha$ and $\beta$ are orthogonal with respect to the weighted inner product
837: $\langle f, g \rangle = \int_{-1}^1 (1+x)^{\alpha} (1-x)^{\beta} f(x) g(x) dx$.
839: Input Parameters:
840: + alpha - the left exponent of the weight
841: . beta - the right exponetn of the weight
842: . npoints - the number of points to evaluate the polynomials at
843: . points - [npoints] array of point coordinates
844: . degree - the maximm degree polynomial space to evaluate, (degree + 1) will be evaluated total.
845: - k - the maximum derivative to evaluate in the jet, (k + 1) will be evaluated total.
847: Output Parameter:
848: . p - an array containing the evaluations of the Jacobi polynomials's jets on the points. the size is (degree + 1) x
849: (k + 1) x npoints, which also describes the order of the dimensions of this three-dimensional array: the first
850: (slowest varying) dimension is polynomial degree; the second dimension is derivative order; the third (fastest
851: varying) dimension is the index of the evaluation point.
853: Level: advanced
855: .seealso: `PetscDTJacobiEval()`, `PetscDTPKDEvalJet()`
856: @*/
857: PetscErrorCode PetscDTJacobiEvalJet(PetscReal alpha, PetscReal beta, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
858: {
859: PetscInt i, j, l;
860: PetscInt *degrees;
861: PetscReal *psingle;
863: PetscFunctionBegin;
864: if (degree == 0) {
865: PetscInt zero = 0;
867: for (i = 0; i <= k; i++) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, 1, &zero, &p[i * npoints]));
868: PetscFunctionReturn(PETSC_SUCCESS);
869: }
870: PetscCall(PetscMalloc1(degree + 1, °rees));
871: PetscCall(PetscMalloc1((degree + 1) * npoints, &psingle));
872: for (i = 0; i <= degree; i++) degrees[i] = i;
873: for (i = 0; i <= k; i++) {
874: PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, degree + 1, degrees, psingle));
875: for (j = 0; j <= degree; j++) {
876: for (l = 0; l < npoints; l++) p[(j * (k + 1) + i) * npoints + l] = psingle[l * (degree + 1) + j];
877: }
878: }
879: PetscCall(PetscFree(psingle));
880: PetscCall(PetscFree(degrees));
881: PetscFunctionReturn(PETSC_SUCCESS);
882: }
884: /*@
885: PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$ at a set of points
886: at points
888: Not Collective
890: Input Parameters:
891: + npoints - number of spatial points to evaluate at
892: . alpha - the left exponent > -1
893: . beta - the right exponent > -1
894: . points - array of locations to evaluate at
895: . ndegree - number of basis degrees to evaluate
896: - degrees - sorted array of degrees to evaluate
898: Output Parameters:
899: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
900: . D - row-oriented derivative evaluation matrix (or NULL)
901: - D2 - row-oriented second derivative evaluation matrix (or NULL)
903: Level: intermediate
905: .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()`
906: @*/
907: PetscErrorCode PetscDTJacobiEval(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *B, PetscReal *D, PetscReal *D2)
908: {
909: PetscFunctionBegin;
910: PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1.");
911: PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1.");
912: if (!npoints || !ndegree) PetscFunctionReturn(PETSC_SUCCESS);
913: if (B) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B));
914: if (D) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D));
915: if (D2) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2));
916: PetscFunctionReturn(PETSC_SUCCESS);
917: }
919: /*@
920: PetscDTLegendreEval - evaluate Legendre polynomials at points
922: Not Collective
924: Input Parameters:
925: + npoints - number of spatial points to evaluate at
926: . points - array of locations to evaluate at
927: . ndegree - number of basis degrees to evaluate
928: - degrees - sorted array of degrees to evaluate
930: Output Parameters:
931: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
932: . D - row-oriented derivative evaluation matrix (or NULL)
933: - D2 - row-oriented second derivative evaluation matrix (or NULL)
935: Level: intermediate
937: .seealso: `PetscDTGaussQuadrature()`
938: @*/
939: PetscErrorCode PetscDTLegendreEval(PetscInt npoints, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *B, PetscReal *D, PetscReal *D2)
940: {
941: PetscFunctionBegin;
942: PetscCall(PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2));
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }
946: /*@
947: PetscDTIndexToGradedOrder - convert an index into a tuple of monomial degrees in a graded order (that is, if the degree sum of tuple x is less than the degree sum of tuple y, then the index of x is smaller than the index of y)
949: Input Parameters:
950: + len - the desired length of the degree tuple
951: - index - the index to convert: should be >= 0
953: Output Parameter:
954: . degtup - will be filled with a tuple of degrees
956: Level: beginner
958: Note:
959: For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
960: acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
961: last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).
963: .seealso: `PetscDTGradedOrderToIndex()`
964: @*/
965: PetscErrorCode PetscDTIndexToGradedOrder(PetscInt len, PetscInt index, PetscInt degtup[])
966: {
967: PetscInt i, total;
968: PetscInt sum;
970: PetscFunctionBeginHot;
971: PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
972: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative");
973: total = 1;
974: sum = 0;
975: while (index >= total) {
976: index -= total;
977: total = (total * (len + sum)) / (sum + 1);
978: sum++;
979: }
980: for (i = 0; i < len; i++) {
981: PetscInt c;
983: degtup[i] = sum;
984: for (c = 0, total = 1; c < sum; c++) {
985: /* going into the loop, total is the number of way to have a tuple of sum exactly c with length len - 1 - i */
986: if (index < total) break;
987: index -= total;
988: total = (total * (len - 1 - i + c)) / (c + 1);
989: degtup[i]--;
990: }
991: sum -= degtup[i];
992: }
993: PetscFunctionReturn(PETSC_SUCCESS);
994: }
996: /*@
997: PetscDTGradedOrderToIndex - convert a tuple into an index in a graded order, the inverse of `PetscDTIndexToGradedOrder()`.
999: Input Parameters:
1000: + len - the length of the degree tuple
1001: - degtup - tuple with this length
1003: Output Parameter:
1004: . index - index in graded order: >= 0
1006: Level: Beginner
1008: Note:
1009: For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
1010: acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
1011: last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).
1013: .seealso: `PetscDTIndexToGradedOrder()`
1014: @*/
1015: PetscErrorCode PetscDTGradedOrderToIndex(PetscInt len, const PetscInt degtup[], PetscInt *index)
1016: {
1017: PetscInt i, idx, sum, total;
1019: PetscFunctionBeginHot;
1020: PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
1021: for (i = 0, sum = 0; i < len; i++) sum += degtup[i];
1022: idx = 0;
1023: total = 1;
1024: for (i = 0; i < sum; i++) {
1025: idx += total;
1026: total = (total * (len + i)) / (i + 1);
1027: }
1028: for (i = 0; i < len - 1; i++) {
1029: PetscInt c;
1031: total = 1;
1032: sum -= degtup[i];
1033: for (c = 0; c < sum; c++) {
1034: idx += total;
1035: total = (total * (len - 1 - i + c)) / (c + 1);
1036: }
1037: }
1038: *index = idx;
1039: PetscFunctionReturn(PETSC_SUCCESS);
1040: }
1042: static PetscBool PKDCite = PETSC_FALSE;
1043: const char PKDCitation[] = "@article{Kirby2010,\n"
1044: " title={Singularity-free evaluation of collapsed-coordinate orthogonal polynomials},\n"
1045: " author={Kirby, Robert C},\n"
1046: " journal={ACM Transactions on Mathematical Software (TOMS)},\n"
1047: " volume={37},\n"
1048: " number={1},\n"
1049: " pages={1--16},\n"
1050: " year={2010},\n"
1051: " publisher={ACM New York, NY, USA}\n}\n";
1053: /*@
1054: PetscDTPKDEvalJet - Evaluate the jet (function and derivatives) of the Proriol-Koornwinder-Dubiner (PKD) basis for
1055: the space of polynomials up to a given degree. The PKD basis is L2-orthonormal on the biunit simplex (which is used
1056: as the reference element for finite elements in PETSc), which makes it a stable basis to use for evaluating
1057: polynomials in that domain.
1059: Input Parameters:
1060: + dim - the number of variables in the multivariate polynomials
1061: . npoints - the number of points to evaluate the polynomials at
1062: . points - [npoints x dim] array of point coordinates
1063: . degree - the degree (sum of degrees on the variables in a monomial) of the polynomial space to evaluate. There are ((dim + degree) choose dim) polynomials in this space.
1064: - k - the maximum order partial derivative to evaluate in the jet. There are (dim + k choose dim) partial derivatives
1065: in the jet. Choosing k = 0 means to evaluate just the function and no derivatives
1067: Output Parameter:
1068: . p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is ((dim + degree)
1069: choose dim) x ((dim + k) choose dim) x npoints, which also describes the order of the dimensions of this
1070: three-dimensional array: the first (slowest varying) dimension is basis function index; the second dimension is jet
1071: index; the third (fastest varying) dimension is the index of the evaluation point.
1073: Level: advanced
1075: Notes:
1076: The ordering of the basis functions, and the ordering of the derivatives in the jet, both follow the graded
1077: ordering of `PetscDTIndexToGradedOrder()` and `PetscDTGradedOrderToIndex()`. For example, in 3D, the polynomial with
1078: leading monomial x^2,y^0,z^1, which has degree tuple (2,0,1), which by `PetscDTGradedOrderToIndex()` has index 12 (it is the 13th basis function in the space);
1079: the partial derivative $\partial_x \partial_z$ has order tuple (1,0,1), appears at index 6 in the jet (it is the 7th partial derivative in the jet).
1081: The implementation uses Kirby's singularity-free evaluation algorithm, https://doi.org/10.1145/1644001.1644006.
1083: .seealso: `PetscDTGradedOrderToIndex()`, `PetscDTIndexToGradedOrder()`, `PetscDTJacobiEvalJet()`
1084: @*/
1085: PetscErrorCode PetscDTPKDEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
1086: {
1087: PetscInt degidx, kidx, d, pt;
1088: PetscInt Nk, Ndeg;
1089: PetscInt *ktup, *degtup;
1090: PetscReal *scales, initscale, scaleexp;
1092: PetscFunctionBegin;
1093: PetscCall(PetscCitationsRegister(PKDCitation, &PKDCite));
1094: PetscCall(PetscDTBinomialInt(dim + k, k, &Nk));
1095: PetscCall(PetscDTBinomialInt(degree + dim, degree, &Ndeg));
1096: PetscCall(PetscMalloc2(dim, °tup, dim, &ktup));
1097: PetscCall(PetscMalloc1(Ndeg, &scales));
1098: initscale = 1.;
1099: if (dim > 1) {
1100: PetscCall(PetscDTBinomial(dim, 2, &scaleexp));
1101: initscale = PetscPowReal(2., scaleexp * 0.5);
1102: }
1103: for (degidx = 0; degidx < Ndeg; degidx++) {
1104: PetscInt e, i;
1105: PetscInt m1idx = -1, m2idx = -1;
1106: PetscInt n;
1107: PetscInt degsum;
1108: PetscReal alpha;
1109: PetscReal cnm1, cnm1x, cnm2;
1110: PetscReal norm;
1112: PetscCall(PetscDTIndexToGradedOrder(dim, degidx, degtup));
1113: for (d = dim - 1; d >= 0; d--)
1114: if (degtup[d]) break;
1115: if (d < 0) { /* constant is 1 everywhere, all derivatives are zero */
1116: scales[degidx] = initscale;
1117: for (e = 0; e < dim; e++) {
1118: PetscCall(PetscDTJacobiNorm(e, 0., 0, &norm));
1119: scales[degidx] /= norm;
1120: }
1121: for (i = 0; i < npoints; i++) p[degidx * Nk * npoints + i] = 1.;
1122: for (i = 0; i < (Nk - 1) * npoints; i++) p[(degidx * Nk + 1) * npoints + i] = 0.;
1123: continue;
1124: }
1125: n = degtup[d];
1126: degtup[d]--;
1127: PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m1idx));
1128: if (degtup[d] > 0) {
1129: degtup[d]--;
1130: PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m2idx));
1131: degtup[d]++;
1132: }
1133: degtup[d]++;
1134: for (e = 0, degsum = 0; e < d; e++) degsum += degtup[e];
1135: alpha = 2 * degsum + d;
1136: PetscDTJacobiRecurrence_Internal(n, alpha, 0., cnm1, cnm1x, cnm2);
1138: scales[degidx] = initscale;
1139: for (e = 0, degsum = 0; e < dim; e++) {
1140: PetscInt f;
1141: PetscReal ealpha;
1142: PetscReal enorm;
1144: ealpha = 2 * degsum + e;
1145: for (f = 0; f < degsum; f++) scales[degidx] *= 2.;
1146: PetscCall(PetscDTJacobiNorm(ealpha, 0., degtup[e], &enorm));
1147: scales[degidx] /= enorm;
1148: degsum += degtup[e];
1149: }
1151: for (pt = 0; pt < npoints; pt++) {
1152: /* compute the multipliers */
1153: PetscReal thetanm1, thetanm1x, thetanm2;
1155: thetanm1x = dim - (d + 1) + 2. * points[pt * dim + d];
1156: for (e = d + 1; e < dim; e++) thetanm1x += points[pt * dim + e];
1157: thetanm1x *= 0.5;
1158: thetanm1 = (2. - (dim - (d + 1)));
1159: for (e = d + 1; e < dim; e++) thetanm1 -= points[pt * dim + e];
1160: thetanm1 *= 0.5;
1161: thetanm2 = thetanm1 * thetanm1;
1163: for (kidx = 0; kidx < Nk; kidx++) {
1164: PetscInt f;
1166: PetscCall(PetscDTIndexToGradedOrder(dim, kidx, ktup));
1167: /* first sum in the same derivative terms */
1168: p[(degidx * Nk + kidx) * npoints + pt] = (cnm1 * thetanm1 + cnm1x * thetanm1x) * p[(m1idx * Nk + kidx) * npoints + pt];
1169: if (m2idx >= 0) p[(degidx * Nk + kidx) * npoints + pt] -= cnm2 * thetanm2 * p[(m2idx * Nk + kidx) * npoints + pt];
1171: for (f = d; f < dim; f++) {
1172: PetscInt km1idx, mplty = ktup[f];
1174: if (!mplty) continue;
1175: ktup[f]--;
1176: PetscCall(PetscDTGradedOrderToIndex(dim, ktup, &km1idx));
1178: /* the derivative of cnm1x * thetanm1x wrt x variable f is 0.5 * cnm1x if f > d otherwise it is cnm1x */
1179: /* the derivative of cnm1 * thetanm1 wrt x variable f is 0 if f == d, otherwise it is -0.5 * cnm1 */
1180: /* the derivative of -cnm2 * thetanm2 wrt x variable f is 0 if f == d, otherwise it is cnm2 * thetanm1 */
1181: if (f > d) {
1182: PetscInt f2;
1184: p[(degidx * Nk + kidx) * npoints + pt] += mplty * 0.5 * (cnm1x - cnm1) * p[(m1idx * Nk + km1idx) * npoints + pt];
1185: if (m2idx >= 0) {
1186: p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm2 * thetanm1 * p[(m2idx * Nk + km1idx) * npoints + pt];
1187: /* second derivatives of -cnm2 * thetanm2 wrt x variable f,f2 is like - 0.5 * cnm2 */
1188: for (f2 = f; f2 < dim; f2++) {
1189: PetscInt km2idx, mplty2 = ktup[f2];
1190: PetscInt factor;
1192: if (!mplty2) continue;
1193: ktup[f2]--;
1194: PetscCall(PetscDTGradedOrderToIndex(dim, ktup, &km2idx));
1196: factor = mplty * mplty2;
1197: if (f == f2) factor /= 2;
1198: p[(degidx * Nk + kidx) * npoints + pt] -= 0.5 * factor * cnm2 * p[(m2idx * Nk + km2idx) * npoints + pt];
1199: ktup[f2]++;
1200: }
1201: }
1202: } else {
1203: p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm1x * p[(m1idx * Nk + km1idx) * npoints + pt];
1204: }
1205: ktup[f]++;
1206: }
1207: }
1208: }
1209: }
1210: for (degidx = 0; degidx < Ndeg; degidx++) {
1211: PetscReal scale = scales[degidx];
1212: PetscInt i;
1214: for (i = 0; i < Nk * npoints; i++) p[degidx * Nk * npoints + i] *= scale;
1215: }
1216: PetscCall(PetscFree(scales));
1217: PetscCall(PetscFree2(degtup, ktup));
1218: PetscFunctionReturn(PETSC_SUCCESS);
1219: }
1221: /*@
1222: PetscDTPTrimmedSize - The size of the trimmed polynomial space of k-forms with a given degree and form degree,
1223: which can be evaluated in `PetscDTPTrimmedEvalJet()`.
1225: Input Parameters:
1226: + dim - the number of variables in the multivariate polynomials
1227: . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space.
1228: - formDegree - the degree of the form
1230: Output Parameter:
1231: - size - The number ((`dim` + `degree`) choose (`dim` + `formDegree`)) x ((`degree` + `formDegree` - 1) choose (`formDegree`))
1233: Level: advanced
1235: .seealso: `PetscDTPTrimmedEvalJet()`
1236: @*/
1237: PetscErrorCode PetscDTPTrimmedSize(PetscInt dim, PetscInt degree, PetscInt formDegree, PetscInt *size)
1238: {
1239: PetscInt Nrk, Nbpt; // number of trimmed polynomials
1241: PetscFunctionBegin;
1242: formDegree = PetscAbsInt(formDegree);
1243: PetscCall(PetscDTBinomialInt(degree + dim, degree + formDegree, &Nbpt));
1244: PetscCall(PetscDTBinomialInt(degree + formDegree - 1, formDegree, &Nrk));
1245: Nbpt *= Nrk;
1246: *size = Nbpt;
1247: PetscFunctionReturn(PETSC_SUCCESS);
1248: }
1250: /* there was a reference implementation based on section 4.4 of Arnold, Falk & Winther (acta numerica, 2006), but it
1251: * was inferior to this implementation */
1252: static PetscErrorCode PetscDTPTrimmedEvalJet_Internal(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[])
1253: {
1254: PetscInt formDegreeOrig = formDegree;
1255: PetscBool formNegative = (formDegreeOrig < 0) ? PETSC_TRUE : PETSC_FALSE;
1257: PetscFunctionBegin;
1258: formDegree = PetscAbsInt(formDegreeOrig);
1259: if (formDegree == 0) {
1260: PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree, jetDegree, p));
1261: PetscFunctionReturn(PETSC_SUCCESS);
1262: }
1263: if (formDegree == dim) {
1264: PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p));
1265: PetscFunctionReturn(PETSC_SUCCESS);
1266: }
1267: PetscInt Nbpt;
1268: PetscCall(PetscDTPTrimmedSize(dim, degree, formDegree, &Nbpt));
1269: PetscInt Nf;
1270: PetscCall(PetscDTBinomialInt(dim, formDegree, &Nf));
1271: PetscInt Nk;
1272: PetscCall(PetscDTBinomialInt(dim + jetDegree, dim, &Nk));
1273: PetscCall(PetscArrayzero(p, Nbpt * Nf * Nk * npoints));
1275: PetscInt Nbpm1; // number of scalar polynomials up to degree - 1;
1276: PetscCall(PetscDTBinomialInt(dim + degree - 1, dim, &Nbpm1));
1277: PetscReal *p_scalar;
1278: PetscCall(PetscMalloc1(Nbpm1 * Nk * npoints, &p_scalar));
1279: PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p_scalar));
1280: PetscInt total = 0;
1281: // First add the full polynomials up to degree - 1 into the basis: take the scalar
1282: // and copy one for each form component
1283: for (PetscInt i = 0; i < Nbpm1; i++) {
1284: const PetscReal *src = &p_scalar[i * Nk * npoints];
1285: for (PetscInt f = 0; f < Nf; f++) {
1286: PetscReal *dest = &p[(total++ * Nf + f) * Nk * npoints];
1287: PetscCall(PetscArraycpy(dest, src, Nk * npoints));
1288: }
1289: }
1290: PetscInt *form_atoms;
1291: PetscCall(PetscMalloc1(formDegree + 1, &form_atoms));
1292: // construct the interior product pattern
1293: PetscInt(*pattern)[3];
1294: PetscInt Nf1; // number of formDegree + 1 forms
1295: PetscCall(PetscDTBinomialInt(dim, formDegree + 1, &Nf1));
1296: PetscInt nnz = Nf1 * (formDegree + 1);
1297: PetscCall(PetscMalloc1(Nf1 * (formDegree + 1), &pattern));
1298: PetscCall(PetscDTAltVInteriorPattern(dim, formDegree + 1, pattern));
1299: PetscReal centroid = (1. - dim) / (dim + 1.);
1300: PetscInt *deriv;
1301: PetscCall(PetscMalloc1(dim, &deriv));
1302: for (PetscInt d = dim; d >= formDegree + 1; d--) {
1303: PetscInt Nfd1; // number of formDegree + 1 forms in dimension d that include dx_0
1304: // (equal to the number of formDegree forms in dimension d-1)
1305: PetscCall(PetscDTBinomialInt(d - 1, formDegree, &Nfd1));
1306: // The number of homogeneous (degree-1) scalar polynomials in d variables
1307: PetscInt Nh;
1308: PetscCall(PetscDTBinomialInt(d - 1 + degree - 1, d - 1, &Nh));
1309: const PetscReal *h_scalar = &p_scalar[(Nbpm1 - Nh) * Nk * npoints];
1310: for (PetscInt b = 0; b < Nh; b++) {
1311: const PetscReal *h_s = &h_scalar[b * Nk * npoints];
1312: for (PetscInt f = 0; f < Nfd1; f++) {
1313: // construct all formDegree+1 forms that start with dx_(dim - d) /\ ...
1314: form_atoms[0] = dim - d;
1315: PetscCall(PetscDTEnumSubset(d - 1, formDegree, f, &form_atoms[1]));
1316: for (PetscInt i = 0; i < formDegree; i++) form_atoms[1 + i] += form_atoms[0] + 1;
1317: PetscInt f_ind; // index of the resulting form
1318: PetscCall(PetscDTSubsetIndex(dim, formDegree + 1, form_atoms, &f_ind));
1319: PetscReal *p_f = &p[total++ * Nf * Nk * npoints];
1320: for (PetscInt nz = 0; nz < nnz; nz++) {
1321: PetscInt i = pattern[nz][0]; // formDegree component
1322: PetscInt j = pattern[nz][1]; // (formDegree + 1) component
1323: PetscInt v = pattern[nz][2]; // coordinate component
1324: PetscReal scale = v < 0 ? -1. : 1.;
1326: i = formNegative ? (Nf - 1 - i) : i;
1327: scale = (formNegative && (i & 1)) ? -scale : scale;
1328: v = v < 0 ? -(v + 1) : v;
1329: if (j != f_ind) continue;
1330: PetscReal *p_i = &p_f[i * Nk * npoints];
1331: for (PetscInt jet = 0; jet < Nk; jet++) {
1332: const PetscReal *h_jet = &h_s[jet * npoints];
1333: PetscReal *p_jet = &p_i[jet * npoints];
1335: for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * h_jet[pt] * (points[pt * dim + v] - centroid);
1336: PetscCall(PetscDTIndexToGradedOrder(dim, jet, deriv));
1337: deriv[v]++;
1338: PetscReal mult = deriv[v];
1339: PetscInt l;
1340: PetscCall(PetscDTGradedOrderToIndex(dim, deriv, &l));
1341: if (l >= Nk) continue;
1342: p_jet = &p_i[l * npoints];
1343: for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * mult * h_jet[pt];
1344: deriv[v]--;
1345: }
1346: }
1347: }
1348: }
1349: }
1350: PetscCheck(total == Nbpt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrectly counted P trimmed polynomials");
1351: PetscCall(PetscFree(deriv));
1352: PetscCall(PetscFree(pattern));
1353: PetscCall(PetscFree(form_atoms));
1354: PetscCall(PetscFree(p_scalar));
1355: PetscFunctionReturn(PETSC_SUCCESS);
1356: }
1358: /*@
1359: PetscDTPTrimmedEvalJet - Evaluate the jet (function and derivatives) of a basis of the trimmed polynomial k-forms up to
1360: a given degree.
1362: Input Parameters:
1363: + dim - the number of variables in the multivariate polynomials
1364: . npoints - the number of points to evaluate the polynomials at
1365: . points - [npoints x dim] array of point coordinates
1366: . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space to evaluate.
1367: There are ((dim + degree) choose (dim + formDegree)) x ((degree + formDegree - 1) choose (formDegree)) polynomials in this space.
1368: (You can use `PetscDTPTrimmedSize()` to compute this size.)
1369: . formDegree - the degree of the form
1370: - jetDegree - the maximum order partial derivative to evaluate in the jet. There are ((dim + jetDegree) choose dim) partial derivatives
1371: in the jet. Choosing jetDegree = 0 means to evaluate just the function and no derivatives
1373: Output Parameter:
1374: . p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is
1375: `PetscDTPTrimmedSize()` x ((dim + formDegree) choose dim) x ((dim + k) choose dim) x npoints,
1376: which also describes the order of the dimensions of this
1377: four-dimensional array:
1378: the first (slowest varying) dimension is basis function index;
1379: the second dimension is component of the form;
1380: the third dimension is jet index;
1381: the fourth (fastest varying) dimension is the index of the evaluation point.
1383: Level: advanced
1385: Notes:
1386: The ordering of the basis functions is not graded, so the basis functions are not nested by degree like `PetscDTPKDEvalJet()`.
1387: The basis functions are not an L2-orthonormal basis on any particular domain.
1389: The implementation is based on the description of the trimmed polynomials up to degree r as
1390: the direct sum of polynomials up to degree (r-1) and the Koszul differential applied to
1391: homogeneous polynomials of degree (r-1).
1393: .seealso: `PetscDTPKDEvalJet()`, `PetscDTPTrimmedSize()`
1394: @*/
1395: PetscErrorCode PetscDTPTrimmedEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[])
1396: {
1397: PetscFunctionBegin;
1398: PetscCall(PetscDTPTrimmedEvalJet_Internal(dim, npoints, points, degree, formDegree, jetDegree, p));
1399: PetscFunctionReturn(PETSC_SUCCESS);
1400: }
1402: /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V
1403: * with lds n; diag and subdiag are overwritten */
1404: static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[], PetscReal eigs[], PetscScalar V[])
1405: {
1406: char jobz = 'V'; /* eigenvalues and eigenvectors */
1407: char range = 'A'; /* all eigenvalues will be found */
1408: PetscReal VL = 0.; /* ignored because range is 'A' */
1409: PetscReal VU = 0.; /* ignored because range is 'A' */
1410: PetscBLASInt IL = 0; /* ignored because range is 'A' */
1411: PetscBLASInt IU = 0; /* ignored because range is 'A' */
1412: PetscReal abstol = 0.; /* unused */
1413: PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */
1414: PetscBLASInt *isuppz;
1415: PetscBLASInt lwork, liwork;
1416: PetscReal workquery;
1417: PetscBLASInt iworkquery;
1418: PetscBLASInt *iwork;
1419: PetscBLASInt info;
1420: PetscReal *work = NULL;
1422: PetscFunctionBegin;
1423: #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1424: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1425: #endif
1426: PetscCall(PetscBLASIntCast(n, &bn));
1427: PetscCall(PetscBLASIntCast(n, &ldz));
1428: #if !defined(PETSC_MISSING_LAPACK_STEGR)
1429: PetscCall(PetscMalloc1(2 * n, &isuppz));
1430: lwork = -1;
1431: liwork = -1;
1432: PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, &workquery, &lwork, &iworkquery, &liwork, &info));
1433: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEGR error");
1434: lwork = (PetscBLASInt)workquery;
1435: liwork = (PetscBLASInt)iworkquery;
1436: PetscCall(PetscMalloc2(lwork, &work, liwork, &iwork));
1437: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1438: PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, work, &lwork, iwork, &liwork, &info));
1439: PetscCall(PetscFPTrapPop());
1440: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEGR error");
1441: PetscCall(PetscFree2(work, iwork));
1442: PetscCall(PetscFree(isuppz));
1443: #elif !defined(PETSC_MISSING_LAPACK_STEQR)
1444: jobz = 'I'; /* Compute eigenvalues and eigenvectors of the
1445: tridiagonal matrix. Z is initialized to the identity
1446: matrix. */
1447: PetscCall(PetscMalloc1(PetscMax(1, 2 * n - 2), &work));
1448: PetscCallBLAS("LAPACKsteqr", LAPACKsteqr_("I", &bn, diag, subdiag, V, &ldz, work, &info));
1449: PetscCall(PetscFPTrapPop());
1450: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEQR error");
1451: PetscCall(PetscFree(work));
1452: PetscCall(PetscArraycpy(eigs, diag, n));
1453: #endif
1454: PetscFunctionReturn(PETSC_SUCCESS);
1455: }
1457: /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi
1458: * quadrature rules on the interval [-1, 1] */
1459: static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw)
1460: {
1461: PetscReal twoab1;
1462: PetscInt m = n - 2;
1463: PetscReal a = alpha + 1.;
1464: PetscReal b = beta + 1.;
1465: PetscReal gra, grb;
1467: PetscFunctionBegin;
1468: twoab1 = PetscPowReal(2., a + b - 1.);
1469: #if defined(PETSC_HAVE_LGAMMA)
1470: grb = PetscExpReal(2. * PetscLGamma(b + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + a + 1.) - (PetscLGamma(m + b + 1) + PetscLGamma(m + a + b + 1.)));
1471: gra = PetscExpReal(2. * PetscLGamma(a + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + b + 1.) - (PetscLGamma(m + a + 1) + PetscLGamma(m + a + b + 1.)));
1472: #else
1473: {
1474: PetscInt alphai = (PetscInt)alpha;
1475: PetscInt betai = (PetscInt)beta;
1477: if ((PetscReal)alphai == alpha && (PetscReal)betai == beta) {
1478: PetscReal binom1, binom2;
1480: PetscCall(PetscDTBinomial(m + b, b, &binom1));
1481: PetscCall(PetscDTBinomial(m + a + b, b, &binom2));
1482: grb = 1. / (binom1 * binom2);
1483: PetscCall(PetscDTBinomial(m + a, a, &binom1));
1484: PetscCall(PetscDTBinomial(m + a + b, a, &binom2));
1485: gra = 1. / (binom1 * binom2);
1486: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
1487: }
1488: #endif
1489: *leftw = twoab1 * grb / b;
1490: *rightw = twoab1 * gra / a;
1491: PetscFunctionReturn(PETSC_SUCCESS);
1492: }
1494: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
1495: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
1496: static inline PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
1497: {
1498: PetscReal pn1, pn2;
1499: PetscReal cnm1, cnm1x, cnm2;
1500: PetscInt k;
1502: PetscFunctionBegin;
1503: if (!n) {
1504: *P = 1.0;
1505: PetscFunctionReturn(PETSC_SUCCESS);
1506: }
1507: PetscDTJacobiRecurrence_Internal(1, a, b, cnm1, cnm1x, cnm2);
1508: pn2 = 1.;
1509: pn1 = cnm1 + cnm1x * x;
1510: if (n == 1) {
1511: *P = pn1;
1512: PetscFunctionReturn(PETSC_SUCCESS);
1513: }
1514: *P = 0.0;
1515: for (k = 2; k < n + 1; ++k) {
1516: PetscDTJacobiRecurrence_Internal(k, a, b, cnm1, cnm1x, cnm2);
1518: *P = (cnm1 + cnm1x * x) * pn1 - cnm2 * pn2;
1519: pn2 = pn1;
1520: pn1 = *P;
1521: }
1522: PetscFunctionReturn(PETSC_SUCCESS);
1523: }
1525: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
1526: static inline PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P)
1527: {
1528: PetscReal nP;
1529: PetscInt i;
1531: PetscFunctionBegin;
1532: *P = 0.0;
1533: if (k > n) PetscFunctionReturn(PETSC_SUCCESS);
1534: PetscCall(PetscDTComputeJacobi(a + k, b + k, n - k, x, &nP));
1535: for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5;
1536: *P = nP;
1537: PetscFunctionReturn(PETSC_SUCCESS);
1538: }
1540: static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[])
1541: {
1542: PetscInt maxIter = 100;
1543: PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON));
1544: PetscReal a1, a6, gf;
1545: PetscInt k;
1547: PetscFunctionBegin;
1549: a1 = PetscPowReal(2.0, a + b + 1);
1550: #if defined(PETSC_HAVE_LGAMMA)
1551: {
1552: PetscReal a2, a3, a4, a5;
1553: a2 = PetscLGamma(a + npoints + 1);
1554: a3 = PetscLGamma(b + npoints + 1);
1555: a4 = PetscLGamma(a + b + npoints + 1);
1556: a5 = PetscLGamma(npoints + 1);
1557: gf = PetscExpReal(a2 + a3 - (a4 + a5));
1558: }
1559: #else
1560: {
1561: PetscInt ia, ib;
1563: ia = (PetscInt)a;
1564: ib = (PetscInt)b;
1565: gf = 1.;
1566: if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */
1567: for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k);
1568: } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */
1569: for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k);
1570: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
1571: }
1572: #endif
1574: a6 = a1 * gf;
1575: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
1576: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
1577: for (k = 0; k < npoints; ++k) {
1578: PetscReal r = PetscCosReal(PETSC_PI * (1. - (4. * k + 3. + 2. * b) / (4. * npoints + 2. * (a + b + 1.)))), dP;
1579: PetscInt j;
1581: if (k > 0) r = 0.5 * (r + x[k - 1]);
1582: for (j = 0; j < maxIter; ++j) {
1583: PetscReal s = 0.0, delta, f, fp;
1584: PetscInt i;
1586: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
1587: PetscCall(PetscDTComputeJacobi(a, b, npoints, r, &f));
1588: PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp));
1589: delta = f / (fp - f * s);
1590: r = r - delta;
1591: if (PetscAbsReal(delta) < eps) break;
1592: }
1593: x[k] = r;
1594: PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP));
1595: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
1596: }
1597: PetscFunctionReturn(PETSC_SUCCESS);
1598: }
1600: /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi
1601: * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */
1602: static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s)
1603: {
1604: PetscInt i;
1606: PetscFunctionBegin;
1607: for (i = 0; i < nPoints; i++) {
1608: PetscReal A, B, C;
1610: PetscDTJacobiRecurrence_Internal(i + 1, a, b, A, B, C);
1611: d[i] = -A / B;
1612: if (i) s[i - 1] *= C / B;
1613: if (i < nPoints - 1) s[i] = 1. / B;
1614: }
1615: PetscFunctionReturn(PETSC_SUCCESS);
1616: }
1618: static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
1619: {
1620: PetscReal mu0;
1621: PetscReal ga, gb, gab;
1622: PetscInt i;
1624: PetscFunctionBegin;
1625: PetscCall(PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite));
1627: #if defined(PETSC_HAVE_TGAMMA)
1628: ga = PetscTGamma(a + 1);
1629: gb = PetscTGamma(b + 1);
1630: gab = PetscTGamma(a + b + 2);
1631: #else
1632: {
1633: PetscInt ia, ib;
1635: ia = (PetscInt)a;
1636: ib = (PetscInt)b;
1637: if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */
1638: PetscCall(PetscDTFactorial(ia, &ga));
1639: PetscCall(PetscDTFactorial(ib, &gb));
1640: PetscCall(PetscDTFactorial(ia + ib + 1, &gab));
1641: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "tgamma() - math routine is unavailable.");
1642: }
1643: #endif
1644: mu0 = PetscPowReal(2., a + b + 1.) * ga * gb / gab;
1646: #if defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1647: {
1648: PetscReal *diag, *subdiag;
1649: PetscScalar *V;
1651: PetscCall(PetscMalloc2(npoints, &diag, npoints, &subdiag));
1652: PetscCall(PetscMalloc1(npoints * npoints, &V));
1653: PetscCall(PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag));
1654: for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]);
1655: PetscCall(PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V));
1656: for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0;
1657: PetscCall(PetscFree(V));
1658: PetscCall(PetscFree2(diag, subdiag));
1659: }
1660: #else
1661: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1662: #endif
1663: { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the
1664: eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that
1665: the eigenvalues are sorted */
1666: PetscBool sorted;
1668: PetscCall(PetscSortedReal(npoints, x, &sorted));
1669: if (!sorted) {
1670: PetscInt *order, i;
1671: PetscReal *tmp;
1673: PetscCall(PetscMalloc2(npoints, &order, npoints, &tmp));
1674: for (i = 0; i < npoints; i++) order[i] = i;
1675: PetscCall(PetscSortRealWithPermutation(npoints, x, order));
1676: PetscCall(PetscArraycpy(tmp, x, npoints));
1677: for (i = 0; i < npoints; i++) x[i] = tmp[order[i]];
1678: PetscCall(PetscArraycpy(tmp, w, npoints));
1679: for (i = 0; i < npoints; i++) w[i] = tmp[order[i]];
1680: PetscCall(PetscFree2(order, tmp));
1681: }
1682: }
1683: PetscFunctionReturn(PETSC_SUCCESS);
1684: }
1686: static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1687: {
1688: PetscFunctionBegin;
1689: PetscCheck(npoints >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive");
1690: /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1691: PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1.");
1692: PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1.");
1694: if (newton) PetscCall(PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w));
1695: else PetscCall(PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w));
1696: if (alpha == beta) { /* symmetrize */
1697: PetscInt i;
1698: for (i = 0; i < (npoints + 1) / 2; i++) {
1699: PetscInt j = npoints - 1 - i;
1700: PetscReal xi = x[i];
1701: PetscReal xj = x[j];
1702: PetscReal wi = w[i];
1703: PetscReal wj = w[j];
1705: x[i] = (xi - xj) / 2.;
1706: x[j] = (xj - xi) / 2.;
1707: w[i] = w[j] = (wi + wj) / 2.;
1708: }
1709: }
1710: PetscFunctionReturn(PETSC_SUCCESS);
1711: }
1713: /*@
1714: PetscDTGaussJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1715: $(x-a)^\alpha (x-b)^\beta$.
1717: Not Collective
1719: Input Parameters:
1720: + npoints - the number of points in the quadrature rule
1721: . a - the left endpoint of the interval
1722: . b - the right endpoint of the interval
1723: . alpha - the left exponent
1724: - beta - the right exponent
1726: Output Parameters:
1727: + x - array of length `npoints`, the locations of the quadrature points
1728: - w - array of length `npoints`, the weights of the quadrature points
1730: Level: intermediate
1732: Note:
1733: This quadrature rule is exact for polynomials up to degree 2*npoints - 1.
1735: .seealso: `PetscDTGaussQuadrature()`
1736: @*/
1737: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1738: {
1739: PetscInt i;
1741: PetscFunctionBegin;
1742: PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal));
1743: if (a != -1. || b != 1.) { /* shift */
1744: for (i = 0; i < npoints; i++) {
1745: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1746: w[i] *= (b - a) / 2.;
1747: }
1748: }
1749: PetscFunctionReturn(PETSC_SUCCESS);
1750: }
1752: static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1753: {
1754: PetscInt i;
1756: PetscFunctionBegin;
1757: PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive");
1758: /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1759: PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1.");
1760: PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1.");
1762: x[0] = -1.;
1763: x[npoints - 1] = 1.;
1764: if (npoints > 2) PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints - 2, alpha + 1., beta + 1., &x[1], &w[1], newton));
1765: for (i = 1; i < npoints - 1; i++) w[i] /= (1. - x[i] * x[i]);
1766: PetscCall(PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints - 1]));
1767: PetscFunctionReturn(PETSC_SUCCESS);
1768: }
1770: /*@
1771: PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1772: $(x-a)^\alpha (x-b)^\beta$, with endpoints a and b included as quadrature points.
1774: Not Collective
1776: Input Parameters:
1777: + npoints - the number of points in the quadrature rule
1778: . a - the left endpoint of the interval
1779: . b - the right endpoint of the interval
1780: . alpha - the left exponent
1781: - beta - the right exponent
1783: Output Parameters:
1784: + x - array of length `npoints`, the locations of the quadrature points
1785: - w - array of length `npoints`, the weights of the quadrature points
1787: Level: intermediate
1789: Note:
1790: This quadrature rule is exact for polynomials up to degree 2*npoints - 3.
1792: .seealso: `PetscDTGaussJacobiQuadrature()`
1793: @*/
1794: PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1795: {
1796: PetscInt i;
1798: PetscFunctionBegin;
1799: PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal));
1800: if (a != -1. || b != 1.) { /* shift */
1801: for (i = 0; i < npoints; i++) {
1802: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1803: w[i] *= (b - a) / 2.;
1804: }
1805: }
1806: PetscFunctionReturn(PETSC_SUCCESS);
1807: }
1809: /*@
1810: PetscDTGaussQuadrature - create Gauss-Legendre quadrature
1812: Not Collective
1814: Input Parameters:
1815: + npoints - number of points
1816: . a - left end of interval (often-1)
1817: - b - right end of interval (often +1)
1819: Output Parameters:
1820: + x - quadrature points
1821: - w - quadrature weights
1823: Level: intermediate
1825: References:
1826: . * - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.
1828: .seealso: `PetscDTLegendreEval()`, `PetscDTGaussJacobiQuadrature()`
1829: @*/
1830: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
1831: {
1832: PetscInt i;
1834: PetscFunctionBegin;
1835: PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal));
1836: if (a != -1. || b != 1.) { /* shift */
1837: for (i = 0; i < npoints; i++) {
1838: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1839: w[i] *= (b - a) / 2.;
1840: }
1841: }
1842: PetscFunctionReturn(PETSC_SUCCESS);
1843: }
1845: /*@C
1846: PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre
1847: nodes of a given size on the domain [-1,1]
1849: Not Collective
1851: Input Parameters:
1852: + n - number of grid nodes
1853: - type - `PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA` or `PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON`
1855: Output Parameters:
1856: + x - quadrature points
1857: - w - quadrature weights
1859: Level: intermediate
1861: Notes:
1862: For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
1863: close enough to the desired solution
1865: These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes
1867: 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
1869: .seealso: `PetscDTGaussQuadrature()`, `PetscGaussLobattoLegendreCreateType`
1871: @*/
1872: PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints, PetscGaussLobattoLegendreCreateType type, PetscReal *x, PetscReal *w)
1873: {
1874: PetscBool newton;
1876: PetscFunctionBegin;
1877: PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide at least 2 grid points per element");
1878: newton = (PetscBool)(type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON);
1879: PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton));
1880: PetscFunctionReturn(PETSC_SUCCESS);
1881: }
1883: /*@
1884: PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
1886: Not Collective
1888: Input Parameters:
1889: + dim - The spatial dimension
1890: . Nc - The number of components
1891: . npoints - number of points in one dimension
1892: . a - left end of interval (often-1)
1893: - b - right end of interval (often +1)
1895: Output Parameter:
1896: . q - A `PetscQuadrature` object
1898: Level: intermediate
1900: .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()`
1901: @*/
1902: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1903: {
1904: DMPolytopeType ct;
1905: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints * PetscSqr(npoints) : PetscSqr(npoints) : npoints;
1906: PetscReal *x, *w, *xw, *ww;
1908: PetscFunctionBegin;
1909: PetscCall(PetscMalloc1(totpoints * dim, &x));
1910: PetscCall(PetscMalloc1(totpoints * Nc, &w));
1911: /* Set up the Golub-Welsch system */
1912: switch (dim) {
1913: case 0:
1914: ct = DM_POLYTOPE_POINT;
1915: PetscCall(PetscFree(x));
1916: PetscCall(PetscFree(w));
1917: PetscCall(PetscMalloc1(1, &x));
1918: PetscCall(PetscMalloc1(Nc, &w));
1919: x[0] = 0.0;
1920: for (PetscInt c = 0; c < Nc; ++c) w[c] = 1.0;
1921: break;
1922: case 1:
1923: ct = DM_POLYTOPE_SEGMENT;
1924: PetscCall(PetscMalloc1(npoints, &ww));
1925: PetscCall(PetscDTGaussQuadrature(npoints, a, b, x, ww));
1926: for (PetscInt i = 0; i < npoints; ++i)
1927: for (PetscInt c = 0; c < Nc; ++c) w[i * Nc + c] = ww[i];
1928: PetscCall(PetscFree(ww));
1929: break;
1930: case 2:
1931: ct = DM_POLYTOPE_QUADRILATERAL;
1932: PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww));
1933: PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww));
1934: for (PetscInt i = 0; i < npoints; ++i) {
1935: for (PetscInt j = 0; j < npoints; ++j) {
1936: x[(i * npoints + j) * dim + 0] = xw[i];
1937: x[(i * npoints + j) * dim + 1] = xw[j];
1938: for (PetscInt c = 0; c < Nc; ++c) w[(i * npoints + j) * Nc + c] = ww[i] * ww[j];
1939: }
1940: }
1941: PetscCall(PetscFree2(xw, ww));
1942: break;
1943: case 3:
1944: ct = DM_POLYTOPE_HEXAHEDRON;
1945: PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww));
1946: PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww));
1947: for (PetscInt i = 0; i < npoints; ++i) {
1948: for (PetscInt j = 0; j < npoints; ++j) {
1949: for (PetscInt k = 0; k < npoints; ++k) {
1950: x[((i * npoints + j) * npoints + k) * dim + 0] = xw[i];
1951: x[((i * npoints + j) * npoints + k) * dim + 1] = xw[j];
1952: x[((i * npoints + j) * npoints + k) * dim + 2] = xw[k];
1953: for (PetscInt c = 0; c < Nc; ++c) w[((i * npoints + j) * npoints + k) * Nc + c] = ww[i] * ww[j] * ww[k];
1954: }
1955: }
1956: }
1957: PetscCall(PetscFree2(xw, ww));
1958: break;
1959: default:
1960: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %" PetscInt_FMT, dim);
1961: }
1962: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
1963: PetscCall(PetscQuadratureSetCellType(*q, ct));
1964: PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1));
1965: PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w));
1966: PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "GaussTensor"));
1967: PetscFunctionReturn(PETSC_SUCCESS);
1968: }
1970: /*@
1971: PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex
1973: Not Collective
1975: Input Parameters:
1976: + dim - The simplex dimension
1977: . Nc - The number of components
1978: . npoints - The number of points in one dimension
1979: . a - left end of interval (often-1)
1980: - b - right end of interval (often +1)
1982: Output Parameter:
1983: . q - A `PetscQuadrature` object
1985: Level: intermediate
1987: Note:
1988: For `dim` == 1, this is Gauss-Legendre quadrature
1990: References:
1991: . * - Karniadakis and Sherwin. FIAT
1993: .seealso: `PetscDTGaussTensorQuadrature()`, `PetscDTGaussQuadrature()`
1994: @*/
1995: PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1996: {
1997: DMPolytopeType ct;
1998: PetscInt totpoints;
1999: PetscReal *p1, *w1;
2000: PetscReal *x, *w;
2002: PetscFunctionBegin;
2003: PetscCheck(!(a != -1.0) && !(b != 1.0), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
2004: switch (dim) {
2005: case 0:
2006: ct = DM_POLYTOPE_POINT;
2007: break;
2008: case 1:
2009: ct = DM_POLYTOPE_SEGMENT;
2010: break;
2011: case 2:
2012: ct = DM_POLYTOPE_TRIANGLE;
2013: break;
2014: case 3:
2015: ct = DM_POLYTOPE_TETRAHEDRON;
2016: break;
2017: default:
2018: ct = DM_POLYTOPE_UNKNOWN;
2019: }
2020: totpoints = 1;
2021: for (PetscInt i = 0; i < dim; ++i) totpoints *= npoints;
2022: PetscCall(PetscMalloc1(totpoints * dim, &x));
2023: PetscCall(PetscMalloc1(totpoints * Nc, &w));
2024: PetscCall(PetscMalloc2(npoints, &p1, npoints, &w1));
2025: for (PetscInt i = 0; i < totpoints * Nc; ++i) w[i] = 1.;
2026: for (PetscInt i = 0, totprev = 1, totrem = totpoints / npoints; i < dim; ++i) {
2027: PetscReal mul;
2029: mul = PetscPowReal(2., -i);
2030: PetscCall(PetscDTGaussJacobiQuadrature(npoints, -1., 1., i, 0.0, p1, w1));
2031: for (PetscInt pt = 0, l = 0; l < totprev; l++) {
2032: for (PetscInt j = 0; j < npoints; j++) {
2033: for (PetscInt m = 0; m < totrem; m++, pt++) {
2034: for (PetscInt k = 0; k < i; k++) x[pt * dim + k] = (x[pt * dim + k] + 1.) * (1. - p1[j]) * 0.5 - 1.;
2035: x[pt * dim + i] = p1[j];
2036: for (PetscInt c = 0; c < Nc; c++) w[pt * Nc + c] *= mul * w1[j];
2037: }
2038: }
2039: }
2040: totprev *= npoints;
2041: totrem /= npoints;
2042: }
2043: PetscCall(PetscFree2(p1, w1));
2044: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
2045: PetscCall(PetscQuadratureSetCellType(*q, ct));
2046: PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1));
2047: PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w));
2048: PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "StroudConical"));
2049: PetscFunctionReturn(PETSC_SUCCESS);
2050: }
2052: static PetscBool MinSymTriQuadCite = PETSC_FALSE;
2053: const char MinSymTriQuadCitation[] = "@article{WitherdenVincent2015,\n"
2054: " title = {On the identification of symmetric quadrature rules for finite element methods},\n"
2055: " journal = {Computers & Mathematics with Applications},\n"
2056: " volume = {69},\n"
2057: " number = {10},\n"
2058: " pages = {1232-1241},\n"
2059: " year = {2015},\n"
2060: " issn = {0898-1221},\n"
2061: " doi = {10.1016/j.camwa.2015.03.017},\n"
2062: " url = {https://www.sciencedirect.com/science/article/pii/S0898122115001224},\n"
2063: " author = {F.D. Witherden and P.E. Vincent},\n"
2064: "}\n";
2066: #include "petscdttriquadrules.h"
2068: static PetscBool MinSymTetQuadCite = PETSC_FALSE;
2069: const char MinSymTetQuadCitation[] = "@article{JaskowiecSukumar2021\n"
2070: " author = {Jaskowiec, Jan and Sukumar, N.},\n"
2071: " title = {High-order symmetric cubature rules for tetrahedra and pyramids},\n"
2072: " journal = {International Journal for Numerical Methods in Engineering},\n"
2073: " volume = {122},\n"
2074: " number = {1},\n"
2075: " pages = {148-171},\n"
2076: " doi = {10.1002/nme.6528},\n"
2077: " url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6528},\n"
2078: " eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6528},\n"
2079: " year = {2021}\n"
2080: "}\n";
2082: #include "petscdttetquadrules.h"
2084: // https://en.wikipedia.org/wiki/Partition_(number_theory)
2085: static PetscErrorCode PetscDTPartitionNumber(PetscInt n, PetscInt *p)
2086: {
2087: // sequence A000041 in the OEIS
2088: const PetscInt partition[] = {1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135, 176, 231, 297, 385, 490, 627, 792, 1002, 1255, 1575, 1958, 2436, 3010, 3718, 4565, 5604};
2089: PetscInt tabulated_max = PETSC_STATIC_ARRAY_LENGTH(partition) - 1;
2091: PetscFunctionBegin;
2092: PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Partition number not defined for negative number %" PetscInt_FMT, n);
2093: // not implementing the pentagonal number recurrence, we don't need partition numbers for n that high
2094: PetscCheck(n <= tabulated_max, PETSC_COMM_SELF, PETSC_ERR_SUP, "Partition numbers only tabulated up to %" PetscInt_FMT ", not computed for %" PetscInt_FMT, tabulated_max, n);
2095: *p = partition[n];
2096: PetscFunctionReturn(PETSC_SUCCESS);
2097: }
2099: /*@
2100: PetscDTSimplexQuadrature - Create a quadrature rule for a simplex that exactly integrates polynomials up to a given degree.
2102: Not Collective
2104: Input Parameters:
2105: + dim - The spatial dimension of the simplex (1 = segment, 2 = triangle, 3 = tetrahedron)
2106: . degree - The largest polynomial degree that is required to be integrated exactly
2107: - type - left end of interval (often-1)
2109: Output Parameter:
2110: . quad - A `PetscQuadrature` object for integration over the biunit simplex
2111: (defined by the bounds $x_i >= -1$ and $\sum_i x_i <= 2 - d$) that is exact for
2112: polynomials up to the given degree
2114: Level: intermediate
2116: .seealso: `PetscDTSimplexQuadratureType`, `PetscDTGaussQuadrature()`, `PetscDTStroudCononicalQuadrature()`, `PetscQuadrature`
2117: @*/
2118: PetscErrorCode PetscDTSimplexQuadrature(PetscInt dim, PetscInt degree, PetscDTSimplexQuadratureType type, PetscQuadrature *quad)
2119: {
2120: PetscDTSimplexQuadratureType orig_type = type;
2122: PetscFunctionBegin;
2123: PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative dimension %" PetscInt_FMT, dim);
2124: PetscCheck(degree >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %" PetscInt_FMT, degree);
2125: if (type == PETSCDTSIMPLEXQUAD_DEFAULT) type = PETSCDTSIMPLEXQUAD_MINSYM;
2126: if (type == PETSCDTSIMPLEXQUAD_CONIC || dim < 2) {
2127: PetscInt points_per_dim = (degree + 2) / 2; // ceil((degree + 1) / 2);
2128: PetscCall(PetscDTStroudConicalQuadrature(dim, 1, points_per_dim, -1, 1, quad));
2129: } else {
2130: DMPolytopeType ct;
2131: PetscInt n = dim + 1;
2132: PetscInt fact = 1;
2133: PetscInt *part, *perm;
2134: PetscInt p = 0;
2135: PetscInt max_degree;
2136: const PetscInt *nodes_per_type = NULL;
2137: const PetscInt *all_num_full_nodes = NULL;
2138: const PetscReal **weights_list = NULL;
2139: const PetscReal **compact_nodes_list = NULL;
2140: const char *citation = NULL;
2141: PetscBool *cited = NULL;
2143: switch (dim) {
2144: case 0:
2145: ct = DM_POLYTOPE_POINT;
2146: break;
2147: case 1:
2148: ct = DM_POLYTOPE_SEGMENT;
2149: break;
2150: case 2:
2151: ct = DM_POLYTOPE_TRIANGLE;
2152: break;
2153: case 3:
2154: ct = DM_POLYTOPE_TETRAHEDRON;
2155: break;
2156: default:
2157: ct = DM_POLYTOPE_UNKNOWN;
2158: }
2159: switch (dim) {
2160: case 2:
2161: cited = &MinSymTriQuadCite;
2162: citation = MinSymTriQuadCitation;
2163: max_degree = PetscDTWVTriQuad_max_degree;
2164: nodes_per_type = PetscDTWVTriQuad_num_orbits;
2165: all_num_full_nodes = PetscDTWVTriQuad_num_nodes;
2166: weights_list = PetscDTWVTriQuad_weights;
2167: compact_nodes_list = PetscDTWVTriQuad_orbits;
2168: break;
2169: case 3:
2170: cited = &MinSymTetQuadCite;
2171: citation = MinSymTetQuadCitation;
2172: max_degree = PetscDTJSTetQuad_max_degree;
2173: nodes_per_type = PetscDTJSTetQuad_num_orbits;
2174: all_num_full_nodes = PetscDTJSTetQuad_num_nodes;
2175: weights_list = PetscDTJSTetQuad_weights;
2176: compact_nodes_list = PetscDTJSTetQuad_orbits;
2177: break;
2178: default:
2179: max_degree = -1;
2180: break;
2181: }
2183: if (degree > max_degree) {
2184: if (orig_type == PETSCDTSIMPLEXQUAD_DEFAULT) {
2185: // fall back to conic
2186: PetscCall(PetscDTSimplexQuadrature(dim, degree, PETSCDTSIMPLEXQUAD_CONIC, quad));
2187: PetscFunctionReturn(PETSC_SUCCESS);
2188: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Minimal symmetric quadrature for dim %" PetscInt_FMT ", degree %" PetscInt_FMT " unsupported", dim, degree);
2189: }
2191: PetscCall(PetscCitationsRegister(citation, cited));
2193: PetscCall(PetscDTPartitionNumber(n, &p));
2194: for (PetscInt d = 2; d <= n; d++) fact *= d;
2196: PetscInt num_full_nodes = all_num_full_nodes[degree];
2197: const PetscReal *all_compact_nodes = compact_nodes_list[degree];
2198: const PetscReal *all_compact_weights = weights_list[degree];
2199: nodes_per_type = &nodes_per_type[p * degree];
2201: PetscReal *points;
2202: PetscReal *counts;
2203: PetscReal *weights;
2204: PetscReal *bary_to_biunit; // row-major transformation of barycentric coordinate to biunit
2205: PetscQuadrature q;
2207: // compute the transformation
2208: PetscCall(PetscMalloc1(n * dim, &bary_to_biunit));
2209: for (PetscInt d = 0; d < dim; d++) {
2210: for (PetscInt b = 0; b < n; b++) bary_to_biunit[d * n + b] = (d == b) ? 1.0 : -1.0;
2211: }
2213: PetscCall(PetscMalloc3(n, &part, n, &perm, n, &counts));
2214: PetscCall(PetscCalloc1(num_full_nodes * dim, &points));
2215: PetscCall(PetscMalloc1(num_full_nodes, &weights));
2217: // (0, 0, ...) is the first partition lexicographically
2218: PetscCall(PetscArrayzero(part, n));
2219: PetscCall(PetscArrayzero(counts, n));
2220: counts[0] = n;
2222: // for each partition
2223: for (PetscInt s = 0, node_offset = 0; s < p; s++) {
2224: PetscInt num_compact_coords = part[n - 1] + 1;
2226: const PetscReal *compact_nodes = all_compact_nodes;
2227: const PetscReal *compact_weights = all_compact_weights;
2228: all_compact_nodes += num_compact_coords * nodes_per_type[s];
2229: all_compact_weights += nodes_per_type[s];
2231: // for every permutation of the vertices
2232: for (PetscInt f = 0; f < fact; f++) {
2233: PetscCall(PetscDTEnumPerm(n, f, perm, NULL));
2235: // check if it is a valid permutation
2236: PetscInt digit;
2237: for (digit = 1; digit < n; digit++) {
2238: // skip permutations that would duplicate a node because it has a smaller symmetry group
2239: if (part[digit - 1] == part[digit] && perm[digit - 1] > perm[digit]) break;
2240: }
2241: if (digit < n) continue;
2243: // create full nodes from this permutation of the compact nodes
2244: PetscReal *full_nodes = &points[node_offset * dim];
2245: PetscReal *full_weights = &weights[node_offset];
2247: PetscCall(PetscArraycpy(full_weights, compact_weights, nodes_per_type[s]));
2248: for (PetscInt b = 0; b < n; b++) {
2249: for (PetscInt d = 0; d < dim; d++) {
2250: for (PetscInt node = 0; node < nodes_per_type[s]; node++) full_nodes[node * dim + d] += bary_to_biunit[d * n + perm[b]] * compact_nodes[node * num_compact_coords + part[b]];
2251: }
2252: }
2253: node_offset += nodes_per_type[s];
2254: }
2256: if (s < p - 1) { // Generate the next partition
2257: /* A partition is described by the number of coordinates that are in
2258: * each set of duplicates (counts) and redundantly by mapping each
2259: * index to its set of duplicates (part)
2260: *
2261: * Counts should always be in nonincreasing order
2262: *
2263: * We want to generate the partitions lexically by part, which means
2264: * finding the last index where count > 1 and reducing by 1.
2265: *
2266: * For the new counts beyond that index, we eagerly assign the remaining
2267: * capacity of the partition to smaller indices (ensures lexical ordering),
2268: * while respecting the nonincreasing invariant of the counts
2269: */
2270: PetscInt last_digit = part[n - 1];
2271: PetscInt last_digit_with_extra = last_digit;
2272: while (counts[last_digit_with_extra] == 1) last_digit_with_extra--;
2273: PetscInt limit = --counts[last_digit_with_extra];
2274: PetscInt total_to_distribute = last_digit - last_digit_with_extra + 1;
2275: for (PetscInt digit = last_digit_with_extra + 1; digit < n; digit++) {
2276: counts[digit] = PetscMin(limit, total_to_distribute);
2277: total_to_distribute -= PetscMin(limit, total_to_distribute);
2278: }
2279: for (PetscInt digit = 0, offset = 0; digit < n; digit++) {
2280: PetscInt count = counts[digit];
2281: for (PetscInt c = 0; c < count; c++) part[offset++] = digit;
2282: }
2283: }
2284: }
2285: PetscCall(PetscFree3(part, perm, counts));
2286: PetscCall(PetscFree(bary_to_biunit));
2287: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &q));
2288: PetscCall(PetscQuadratureSetCellType(q, ct));
2289: PetscCall(PetscQuadratureSetOrder(q, degree));
2290: PetscCall(PetscQuadratureSetData(q, dim, 1, num_full_nodes, points, weights));
2291: *quad = q;
2292: }
2293: PetscFunctionReturn(PETSC_SUCCESS);
2294: }
2296: /*@
2297: PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell
2299: Not Collective
2301: Input Parameters:
2302: + dim - The cell dimension
2303: . level - The number of points in one dimension, 2^l
2304: . a - left end of interval (often-1)
2305: - b - right end of interval (often +1)
2307: Output Parameter:
2308: . q - A `PetscQuadrature` object
2310: Level: intermediate
2312: .seealso: `PetscDTGaussTensorQuadrature()`, `PetscQuadrature`
2313: @*/
2314: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
2315: {
2316: DMPolytopeType ct;
2317: const PetscInt p = 16; /* Digits of precision in the evaluation */
2318: const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */
2319: const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */
2320: const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */
2321: PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */
2322: PetscReal wk = 0.5 * PETSC_PI; /* Quadrature weight at x_k */
2323: PetscReal *x, *w;
2324: PetscInt K, k, npoints;
2326: PetscFunctionBegin;
2327: PetscCheck(dim <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not yet implemented", dim);
2328: PetscCheck(level, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
2329: switch (dim) {
2330: case 0:
2331: ct = DM_POLYTOPE_POINT;
2332: break;
2333: case 1:
2334: ct = DM_POLYTOPE_SEGMENT;
2335: break;
2336: case 2:
2337: ct = DM_POLYTOPE_QUADRILATERAL;
2338: break;
2339: case 3:
2340: ct = DM_POLYTOPE_HEXAHEDRON;
2341: break;
2342: default:
2343: ct = DM_POLYTOPE_UNKNOWN;
2344: }
2345: /* Find K such that the weights are < 32 digits of precision */
2346: for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2 * p; ++K) wk = 0.5 * h * PETSC_PI * PetscCoshReal(K * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(K * h)));
2347: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
2348: PetscCall(PetscQuadratureSetCellType(*q, ct));
2349: PetscCall(PetscQuadratureSetOrder(*q, 2 * K + 1));
2350: npoints = 2 * K - 1;
2351: PetscCall(PetscMalloc1(npoints * dim, &x));
2352: PetscCall(PetscMalloc1(npoints, &w));
2353: /* Center term */
2354: x[0] = beta;
2355: w[0] = 0.5 * alpha * PETSC_PI;
2356: for (k = 1; k < K; ++k) {
2357: wk = 0.5 * alpha * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2358: xk = PetscTanhReal(0.5 * PETSC_PI * PetscSinhReal(k * h));
2359: x[2 * k - 1] = -alpha * xk + beta;
2360: w[2 * k - 1] = wk;
2361: x[2 * k + 0] = alpha * xk + beta;
2362: w[2 * k + 0] = wk;
2363: }
2364: PetscCall(PetscQuadratureSetData(*q, dim, 1, npoints, x, w));
2365: PetscFunctionReturn(PETSC_SUCCESS);
2366: }
2368: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2369: {
2370: const PetscInt p = 16; /* Digits of precision in the evaluation */
2371: const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */
2372: const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */
2373: PetscReal h = 1.0; /* Step size, length between x_k */
2374: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
2375: PetscReal osum = 0.0; /* Integral on last level */
2376: PetscReal psum = 0.0; /* Integral on the level before the last level */
2377: PetscReal sum; /* Integral on current level */
2378: PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
2379: PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
2380: PetscReal wk; /* Quadrature weight at x_k */
2381: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
2382: PetscInt d; /* Digits of precision in the integral */
2384: PetscFunctionBegin;
2385: PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
2386: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2387: /* Center term */
2388: func(&beta, ctx, &lval);
2389: sum = 0.5 * alpha * PETSC_PI * lval;
2390: /* */
2391: do {
2392: PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
2393: PetscInt k = 1;
2395: ++l;
2396: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */
2397: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
2398: psum = osum;
2399: osum = sum;
2400: h *= 0.5;
2401: sum *= 0.5;
2402: do {
2403: wk = 0.5 * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2404: yk = 1.0 / (PetscExpReal(0.5 * PETSC_PI * PetscSinhReal(k * h)) * PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2405: lx = -alpha * (1.0 - yk) + beta;
2406: rx = alpha * (1.0 - yk) + beta;
2407: func(&lx, ctx, &lval);
2408: func(&rx, ctx, &rval);
2409: lterm = alpha * wk * lval;
2410: maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
2411: sum += lterm;
2412: rterm = alpha * wk * rval;
2413: maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
2414: sum += rterm;
2415: ++k;
2416: /* Only need to evaluate every other point on refined levels */
2417: if (l != 1) ++k;
2418: } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */
2420: d1 = PetscLog10Real(PetscAbsReal(sum - osum));
2421: d2 = PetscLog10Real(PetscAbsReal(sum - psum));
2422: d3 = PetscLog10Real(maxTerm) - p;
2423: if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
2424: else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
2425: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4)));
2426: } while (d < digits && l < 12);
2427: *sol = sum;
2428: PetscCall(PetscFPTrapPop());
2429: PetscFunctionReturn(PETSC_SUCCESS);
2430: }
2432: #if defined(PETSC_HAVE_MPFR)
2433: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2434: {
2435: const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */
2436: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
2437: mpfr_t alpha; /* Half-width of the integration interval */
2438: mpfr_t beta; /* Center of the integration interval */
2439: mpfr_t h; /* Step size, length between x_k */
2440: mpfr_t osum; /* Integral on last level */
2441: mpfr_t psum; /* Integral on the level before the last level */
2442: mpfr_t sum; /* Integral on current level */
2443: mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
2444: mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
2445: mpfr_t wk; /* Quadrature weight at x_k */
2446: PetscReal lval, rval, rtmp; /* Terms in the quadature sum to the left and right of 0 */
2447: PetscInt d; /* Digits of precision in the integral */
2448: mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;
2450: PetscFunctionBegin;
2451: PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
2452: /* Create high precision storage */
2453: 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);
2454: /* Initialization */
2455: mpfr_set_d(alpha, 0.5 * (b - a), MPFR_RNDN);
2456: mpfr_set_d(beta, 0.5 * (b + a), MPFR_RNDN);
2457: mpfr_set_d(osum, 0.0, MPFR_RNDN);
2458: mpfr_set_d(psum, 0.0, MPFR_RNDN);
2459: mpfr_set_d(h, 1.0, MPFR_RNDN);
2460: mpfr_const_pi(pi2, MPFR_RNDN);
2461: mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
2462: /* Center term */
2463: rtmp = 0.5 * (b + a);
2464: func(&rtmp, ctx, &lval);
2465: mpfr_set(sum, pi2, MPFR_RNDN);
2466: mpfr_mul(sum, sum, alpha, MPFR_RNDN);
2467: mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
2468: /* */
2469: do {
2470: PetscReal d1, d2, d3, d4;
2471: PetscInt k = 1;
2473: ++l;
2474: mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
2475: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */
2476: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
2477: mpfr_set(psum, osum, MPFR_RNDN);
2478: mpfr_set(osum, sum, MPFR_RNDN);
2479: mpfr_mul_d(h, h, 0.5, MPFR_RNDN);
2480: mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
2481: do {
2482: mpfr_set_si(kh, k, MPFR_RNDN);
2483: mpfr_mul(kh, kh, h, MPFR_RNDN);
2484: /* Weight */
2485: mpfr_set(wk, h, MPFR_RNDN);
2486: mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
2487: mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
2488: mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
2489: mpfr_cosh(tmp, msinh, MPFR_RNDN);
2490: mpfr_sqr(tmp, tmp, MPFR_RNDN);
2491: mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
2492: mpfr_div(wk, wk, tmp, MPFR_RNDN);
2493: /* Abscissa */
2494: mpfr_set_d(yk, 1.0, MPFR_RNDZ);
2495: mpfr_cosh(tmp, msinh, MPFR_RNDN);
2496: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
2497: mpfr_exp(tmp, msinh, MPFR_RNDN);
2498: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
2499: /* Quadrature points */
2500: mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
2501: mpfr_mul(lx, lx, alpha, MPFR_RNDU);
2502: mpfr_add(lx, lx, beta, MPFR_RNDU);
2503: mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
2504: mpfr_mul(rx, rx, alpha, MPFR_RNDD);
2505: mpfr_add(rx, rx, beta, MPFR_RNDD);
2506: /* Evaluation */
2507: rtmp = mpfr_get_d(lx, MPFR_RNDU);
2508: func(&rtmp, ctx, &lval);
2509: rtmp = mpfr_get_d(rx, MPFR_RNDD);
2510: func(&rtmp, ctx, &rval);
2511: /* Update */
2512: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
2513: mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
2514: mpfr_add(sum, sum, tmp, MPFR_RNDN);
2515: mpfr_abs(tmp, tmp, MPFR_RNDN);
2516: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
2517: mpfr_set(curTerm, tmp, MPFR_RNDN);
2518: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
2519: mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
2520: mpfr_add(sum, sum, tmp, MPFR_RNDN);
2521: mpfr_abs(tmp, tmp, MPFR_RNDN);
2522: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
2523: mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
2524: ++k;
2525: /* Only need to evaluate every other point on refined levels */
2526: if (l != 1) ++k;
2527: mpfr_log10(tmp, wk, MPFR_RNDN);
2528: mpfr_abs(tmp, tmp, MPFR_RNDN);
2529: } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor * digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
2530: mpfr_sub(tmp, sum, osum, MPFR_RNDN);
2531: mpfr_abs(tmp, tmp, MPFR_RNDN);
2532: mpfr_log10(tmp, tmp, MPFR_RNDN);
2533: d1 = mpfr_get_d(tmp, MPFR_RNDN);
2534: mpfr_sub(tmp, sum, psum, MPFR_RNDN);
2535: mpfr_abs(tmp, tmp, MPFR_RNDN);
2536: mpfr_log10(tmp, tmp, MPFR_RNDN);
2537: d2 = mpfr_get_d(tmp, MPFR_RNDN);
2538: mpfr_log10(tmp, maxTerm, MPFR_RNDN);
2539: d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
2540: mpfr_log10(tmp, curTerm, MPFR_RNDN);
2541: d4 = mpfr_get_d(tmp, MPFR_RNDN);
2542: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4)));
2543: } while (d < digits && l < 8);
2544: *sol = mpfr_get_d(sum, MPFR_RNDN);
2545: /* Cleanup */
2546: mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
2547: PetscFunctionReturn(PETSC_SUCCESS);
2548: }
2549: #else
2551: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2552: {
2553: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
2554: }
2555: #endif
2557: /*@
2558: PetscDTTensorQuadratureCreate - create the tensor product quadrature from two lower-dimensional quadratures
2560: Not Collective
2562: Input Parameters:
2563: + q1 - The first quadrature
2564: - q2 - The second quadrature
2566: Output Parameter:
2567: . q - A `PetscQuadrature` object
2569: Level: intermediate
2571: .seealso: `PetscQuadrature`, `PetscDTGaussTensorQuadrature()`
2572: @*/
2573: PetscErrorCode PetscDTTensorQuadratureCreate(PetscQuadrature q1, PetscQuadrature q2, PetscQuadrature *q)
2574: {
2575: DMPolytopeType ct1, ct2, ct;
2576: const PetscReal *x1, *w1, *x2, *w2;
2577: PetscReal *x, *w;
2578: PetscInt dim1, Nc1, Np1, order1, qa, d1;
2579: PetscInt dim2, Nc2, Np2, order2, qb, d2;
2580: PetscInt dim, Nc, Np, order, qc, d;
2582: PetscFunctionBegin;
2586: PetscCall(PetscQuadratureGetOrder(q1, &order1));
2587: PetscCall(PetscQuadratureGetOrder(q2, &order2));
2588: PetscCheck(order1 == order2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Order1 %" PetscInt_FMT " != %" PetscInt_FMT " Order2", order1, order2);
2589: PetscCall(PetscQuadratureGetData(q1, &dim1, &Nc1, &Np1, &x1, &w1));
2590: PetscCall(PetscQuadratureGetCellType(q1, &ct1));
2591: PetscCall(PetscQuadratureGetData(q2, &dim2, &Nc2, &Np2, &x2, &w2));
2592: PetscCall(PetscQuadratureGetCellType(q2, &ct2));
2593: PetscCheck(Nc1 == Nc2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "NumComp1 %" PetscInt_FMT " != %" PetscInt_FMT " NumComp2", Nc1, Nc2);
2595: switch (ct1) {
2596: case DM_POLYTOPE_POINT:
2597: ct = ct2;
2598: break;
2599: case DM_POLYTOPE_SEGMENT:
2600: switch (ct2) {
2601: case DM_POLYTOPE_POINT:
2602: ct = ct1;
2603: break;
2604: case DM_POLYTOPE_SEGMENT:
2605: ct = DM_POLYTOPE_QUADRILATERAL;
2606: break;
2607: case DM_POLYTOPE_TRIANGLE:
2608: ct = DM_POLYTOPE_TRI_PRISM;
2609: break;
2610: case DM_POLYTOPE_QUADRILATERAL:
2611: ct = DM_POLYTOPE_HEXAHEDRON;
2612: break;
2613: case DM_POLYTOPE_TETRAHEDRON:
2614: ct = DM_POLYTOPE_UNKNOWN;
2615: break;
2616: case DM_POLYTOPE_HEXAHEDRON:
2617: ct = DM_POLYTOPE_UNKNOWN;
2618: break;
2619: default:
2620: ct = DM_POLYTOPE_UNKNOWN;
2621: }
2622: break;
2623: case DM_POLYTOPE_TRIANGLE:
2624: switch (ct2) {
2625: case DM_POLYTOPE_POINT:
2626: ct = ct1;
2627: break;
2628: case DM_POLYTOPE_SEGMENT:
2629: ct = DM_POLYTOPE_TRI_PRISM;
2630: break;
2631: case DM_POLYTOPE_TRIANGLE:
2632: ct = DM_POLYTOPE_UNKNOWN;
2633: break;
2634: case DM_POLYTOPE_QUADRILATERAL:
2635: ct = DM_POLYTOPE_UNKNOWN;
2636: break;
2637: case DM_POLYTOPE_TETRAHEDRON:
2638: ct = DM_POLYTOPE_UNKNOWN;
2639: break;
2640: case DM_POLYTOPE_HEXAHEDRON:
2641: ct = DM_POLYTOPE_UNKNOWN;
2642: break;
2643: default:
2644: ct = DM_POLYTOPE_UNKNOWN;
2645: }
2646: break;
2647: case DM_POLYTOPE_QUADRILATERAL:
2648: switch (ct2) {
2649: case DM_POLYTOPE_POINT:
2650: ct = ct1;
2651: break;
2652: case DM_POLYTOPE_SEGMENT:
2653: ct = DM_POLYTOPE_HEXAHEDRON;
2654: break;
2655: case DM_POLYTOPE_TRIANGLE:
2656: ct = DM_POLYTOPE_UNKNOWN;
2657: break;
2658: case DM_POLYTOPE_QUADRILATERAL:
2659: ct = DM_POLYTOPE_UNKNOWN;
2660: break;
2661: case DM_POLYTOPE_TETRAHEDRON:
2662: ct = DM_POLYTOPE_UNKNOWN;
2663: break;
2664: case DM_POLYTOPE_HEXAHEDRON:
2665: ct = DM_POLYTOPE_UNKNOWN;
2666: break;
2667: default:
2668: ct = DM_POLYTOPE_UNKNOWN;
2669: }
2670: break;
2671: case DM_POLYTOPE_TETRAHEDRON:
2672: switch (ct2) {
2673: case DM_POLYTOPE_POINT:
2674: ct = ct1;
2675: break;
2676: case DM_POLYTOPE_SEGMENT:
2677: ct = DM_POLYTOPE_UNKNOWN;
2678: break;
2679: case DM_POLYTOPE_TRIANGLE:
2680: ct = DM_POLYTOPE_UNKNOWN;
2681: break;
2682: case DM_POLYTOPE_QUADRILATERAL:
2683: ct = DM_POLYTOPE_UNKNOWN;
2684: break;
2685: case DM_POLYTOPE_TETRAHEDRON:
2686: ct = DM_POLYTOPE_UNKNOWN;
2687: break;
2688: case DM_POLYTOPE_HEXAHEDRON:
2689: ct = DM_POLYTOPE_UNKNOWN;
2690: break;
2691: default:
2692: ct = DM_POLYTOPE_UNKNOWN;
2693: }
2694: break;
2695: case DM_POLYTOPE_HEXAHEDRON:
2696: switch (ct2) {
2697: case DM_POLYTOPE_POINT:
2698: ct = ct1;
2699: break;
2700: case DM_POLYTOPE_SEGMENT:
2701: ct = DM_POLYTOPE_UNKNOWN;
2702: break;
2703: case DM_POLYTOPE_TRIANGLE:
2704: ct = DM_POLYTOPE_UNKNOWN;
2705: break;
2706: case DM_POLYTOPE_QUADRILATERAL:
2707: ct = DM_POLYTOPE_UNKNOWN;
2708: break;
2709: case DM_POLYTOPE_TETRAHEDRON:
2710: ct = DM_POLYTOPE_UNKNOWN;
2711: break;
2712: case DM_POLYTOPE_HEXAHEDRON:
2713: ct = DM_POLYTOPE_UNKNOWN;
2714: break;
2715: default:
2716: ct = DM_POLYTOPE_UNKNOWN;
2717: }
2718: break;
2719: default:
2720: ct = DM_POLYTOPE_UNKNOWN;
2721: }
2722: dim = dim1 + dim2;
2723: Nc = Nc1;
2724: Np = Np1 * Np2;
2725: order = order1;
2726: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
2727: PetscCall(PetscQuadratureSetCellType(*q, ct));
2728: PetscCall(PetscQuadratureSetOrder(*q, order));
2729: PetscCall(PetscMalloc1(Np * dim, &x));
2730: PetscCall(PetscMalloc1(Np, &w));
2731: for (qa = 0, qc = 0; qa < Np1; ++qa) {
2732: for (qb = 0; qb < Np2; ++qb, ++qc) {
2733: for (d1 = 0, d = 0; d1 < dim1; ++d1, ++d) x[qc * dim + d] = x1[qa * dim1 + d1];
2734: for (d2 = 0; d2 < dim2; ++d2, ++d) x[qc * dim + d] = x2[qb * dim2 + d2];
2735: w[qc] = w1[qa] * w2[qb];
2736: }
2737: }
2738: PetscCall(PetscQuadratureSetData(*q, dim, Nc, Np, x, w));
2739: PetscFunctionReturn(PETSC_SUCCESS);
2740: }
2742: /* Overwrites A. Can only handle full-rank problems with m>=n
2743: A in column-major format
2744: Ainv in row-major format
2745: tau has length m
2746: worksize must be >= max(1,n)
2747: */
2748: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m, PetscInt mstride, PetscInt n, PetscReal *A_in, PetscReal *Ainv_out, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2749: {
2750: PetscBLASInt M, N, K, lda, ldb, ldwork, info;
2751: PetscScalar *A, *Ainv, *R, *Q, Alpha;
2753: PetscFunctionBegin;
2754: #if defined(PETSC_USE_COMPLEX)
2755: {
2756: PetscInt i, j;
2757: PetscCall(PetscMalloc2(m * n, &A, m * n, &Ainv));
2758: for (j = 0; j < n; j++) {
2759: for (i = 0; i < m; i++) A[i + m * j] = A_in[i + mstride * j];
2760: }
2761: mstride = m;
2762: }
2763: #else
2764: A = A_in;
2765: Ainv = Ainv_out;
2766: #endif
2768: PetscCall(PetscBLASIntCast(m, &M));
2769: PetscCall(PetscBLASIntCast(n, &N));
2770: PetscCall(PetscBLASIntCast(mstride, &lda));
2771: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2772: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2773: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
2774: PetscCall(PetscFPTrapPop());
2775: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
2776: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2778: /* Extract an explicit representation of Q */
2779: Q = Ainv;
2780: PetscCall(PetscArraycpy(Q, A, mstride * n));
2781: K = N; /* full rank */
2782: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
2783: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");
2785: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2786: Alpha = 1.0;
2787: ldb = lda;
2788: PetscCallBLAS("BLAStrsm", BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb));
2789: /* Ainv is Q, overwritten with inverse */
2791: #if defined(PETSC_USE_COMPLEX)
2792: {
2793: PetscInt i;
2794: for (i = 0; i < m * n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
2795: PetscCall(PetscFree2(A, Ainv));
2796: }
2797: #endif
2798: PetscFunctionReturn(PETSC_SUCCESS);
2799: }
2801: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
2802: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval, const PetscReal *x, PetscInt ndegree, const PetscInt *degrees, PetscBool Transpose, PetscReal *B)
2803: {
2804: PetscReal *Bv;
2805: PetscInt i, j;
2807: PetscFunctionBegin;
2808: PetscCall(PetscMalloc1((ninterval + 1) * ndegree, &Bv));
2809: /* Point evaluation of L_p on all the source vertices */
2810: PetscCall(PetscDTLegendreEval(ninterval + 1, x, ndegree, degrees, Bv, NULL, NULL));
2811: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
2812: for (i = 0; i < ninterval; i++) {
2813: for (j = 0; j < ndegree; j++) {
2814: if (Transpose) B[i + ninterval * j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j];
2815: else B[i * ndegree + j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j];
2816: }
2817: }
2818: PetscCall(PetscFree(Bv));
2819: PetscFunctionReturn(PETSC_SUCCESS);
2820: }
2822: /*@
2823: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
2825: Not Collective
2827: Input Parameters:
2828: + degree - degree of reconstruction polynomial
2829: . nsource - number of source intervals
2830: . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
2831: . ntarget - number of target intervals
2832: - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
2834: Output Parameter:
2835: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
2837: Level: advanced
2839: .seealso: `PetscDTLegendreEval()`
2840: @*/
2841: PetscErrorCode PetscDTReconstructPoly(PetscInt degree, PetscInt nsource, const PetscReal *sourcex, PetscInt ntarget, const PetscReal *targetx, PetscReal *R)
2842: {
2843: PetscInt i, j, k, *bdegrees, worksize;
2844: PetscReal xmin, xmax, center, hscale, *sourcey, *targety, *Bsource, *Bsinv, *Btarget;
2845: PetscScalar *tau, *work;
2847: PetscFunctionBegin;
2851: PetscCheck(degree < nsource, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Reconstruction degree %" PetscInt_FMT " must be less than number of source intervals %" PetscInt_FMT, degree, nsource);
2852: if (PetscDefined(USE_DEBUG)) {
2853: for (i = 0; i < nsource; i++) PetscCheck(sourcex[i] < sourcex[i + 1], PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Source interval %" PetscInt_FMT " has negative orientation (%g,%g)", i, (double)sourcex[i], (double)sourcex[i + 1]);
2854: for (i = 0; i < ntarget; i++) PetscCheck(targetx[i] < targetx[i + 1], PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Target interval %" PetscInt_FMT " has negative orientation (%g,%g)", i, (double)targetx[i], (double)targetx[i + 1]);
2855: }
2856: xmin = PetscMin(sourcex[0], targetx[0]);
2857: xmax = PetscMax(sourcex[nsource], targetx[ntarget]);
2858: center = (xmin + xmax) / 2;
2859: hscale = (xmax - xmin) / 2;
2860: worksize = nsource;
2861: PetscCall(PetscMalloc4(degree + 1, &bdegrees, nsource + 1, &sourcey, nsource * (degree + 1), &Bsource, worksize, &work));
2862: PetscCall(PetscMalloc4(nsource, &tau, nsource * (degree + 1), &Bsinv, ntarget + 1, &targety, ntarget * (degree + 1), &Btarget));
2863: for (i = 0; i <= nsource; i++) sourcey[i] = (sourcex[i] - center) / hscale;
2864: for (i = 0; i <= degree; i++) bdegrees[i] = i + 1;
2865: PetscCall(PetscDTLegendreIntegrate(nsource, sourcey, degree + 1, bdegrees, PETSC_TRUE, Bsource));
2866: PetscCall(PetscDTPseudoInverseQR(nsource, nsource, degree + 1, Bsource, Bsinv, tau, nsource, work));
2867: for (i = 0; i <= ntarget; i++) targety[i] = (targetx[i] - center) / hscale;
2868: PetscCall(PetscDTLegendreIntegrate(ntarget, targety, degree + 1, bdegrees, PETSC_FALSE, Btarget));
2869: for (i = 0; i < ntarget; i++) {
2870: PetscReal rowsum = 0;
2871: for (j = 0; j < nsource; j++) {
2872: PetscReal sum = 0;
2873: for (k = 0; k < degree + 1; k++) sum += Btarget[i * (degree + 1) + k] * Bsinv[k * nsource + j];
2874: R[i * nsource + j] = sum;
2875: rowsum += sum;
2876: }
2877: for (j = 0; j < nsource; j++) R[i * nsource + j] /= rowsum; /* normalize each row */
2878: }
2879: PetscCall(PetscFree4(bdegrees, sourcey, Bsource, work));
2880: PetscCall(PetscFree4(tau, Bsinv, targety, Btarget));
2881: PetscFunctionReturn(PETSC_SUCCESS);
2882: }
2884: /*@C
2885: PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points
2887: Not Collective
2889: Input Parameters:
2890: + n - the number of GLL nodes
2891: . nodes - the GLL nodes
2892: . weights - the GLL weights
2893: - f - the function values at the nodes
2895: Output Parameter:
2896: . in - the value of the integral
2898: Level: beginner
2900: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`
2901: @*/
2902: PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n, PetscReal *nodes, PetscReal *weights, const PetscReal *f, PetscReal *in)
2903: {
2904: PetscInt i;
2906: PetscFunctionBegin;
2907: *in = 0.;
2908: for (i = 0; i < n; i++) *in += f[i] * f[i] * weights[i];
2909: PetscFunctionReturn(PETSC_SUCCESS);
2910: }
2912: /*@C
2913: PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element
2915: Not Collective
2917: Input Parameters:
2918: + n - the number of GLL nodes
2919: . nodes - the GLL nodes
2920: - weights - the GLL weights
2922: Output Parameter:
2923: . A - the stiffness element
2925: Level: beginner
2927: Notes:
2928: Destroy this with `PetscGaussLobattoLegendreElementLaplacianDestroy()`
2930: 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)
2932: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`
2933: @*/
2934: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
2935: {
2936: PetscReal **A;
2937: const PetscReal *gllnodes = nodes;
2938: const PetscInt p = n - 1;
2939: PetscReal z0, z1, z2 = -1, x, Lpj, Lpr;
2940: PetscInt i, j, nn, r;
2942: PetscFunctionBegin;
2943: PetscCall(PetscMalloc1(n, &A));
2944: PetscCall(PetscMalloc1(n * n, &A[0]));
2945: for (i = 1; i < n; i++) A[i] = A[i - 1] + n;
2947: for (j = 1; j < p; j++) {
2948: x = gllnodes[j];
2949: z0 = 1.;
2950: z1 = x;
2951: for (nn = 1; nn < p; nn++) {
2952: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2953: z0 = z1;
2954: z1 = z2;
2955: }
2956: Lpj = z2;
2957: for (r = 1; r < p; r++) {
2958: if (r == j) {
2959: A[j][j] = 2. / (3. * (1. - gllnodes[j] * gllnodes[j]) * Lpj * Lpj);
2960: } else {
2961: x = gllnodes[r];
2962: z0 = 1.;
2963: z1 = x;
2964: for (nn = 1; nn < p; nn++) {
2965: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2966: z0 = z1;
2967: z1 = z2;
2968: }
2969: Lpr = z2;
2970: A[r][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * Lpr * (gllnodes[j] - gllnodes[r]) * (gllnodes[j] - gllnodes[r]));
2971: }
2972: }
2973: }
2974: for (j = 1; j < p + 1; j++) {
2975: x = gllnodes[j];
2976: z0 = 1.;
2977: z1 = x;
2978: for (nn = 1; nn < p; nn++) {
2979: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2980: z0 = z1;
2981: z1 = z2;
2982: }
2983: Lpj = z2;
2984: A[j][0] = 4. * PetscPowRealInt(-1., p) / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. + gllnodes[j]) * (1. + gllnodes[j]));
2985: A[0][j] = A[j][0];
2986: }
2987: for (j = 0; j < p; j++) {
2988: x = gllnodes[j];
2989: z0 = 1.;
2990: z1 = x;
2991: for (nn = 1; nn < p; nn++) {
2992: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2993: z0 = z1;
2994: z1 = z2;
2995: }
2996: Lpj = z2;
2998: A[p][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. - gllnodes[j]) * (1. - gllnodes[j]));
2999: A[j][p] = A[p][j];
3000: }
3001: A[0][0] = 0.5 + (((PetscReal)p) * (((PetscReal)p) + 1.) - 2.) / 6.;
3002: A[p][p] = A[0][0];
3003: *AA = A;
3004: PetscFunctionReturn(PETSC_SUCCESS);
3005: }
3007: /*@C
3008: PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element created with `PetscGaussLobattoLegendreElementLaplacianCreate()`
3010: Not Collective
3012: Input Parameters:
3013: + n - the number of GLL nodes
3014: . nodes - the GLL nodes
3015: . weights - the GLL weightss
3016: - A - the stiffness element
3018: Level: beginner
3020: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`
3021: @*/
3022: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3023: {
3024: PetscFunctionBegin;
3025: PetscCall(PetscFree((*AA)[0]));
3026: PetscCall(PetscFree(*AA));
3027: *AA = NULL;
3028: PetscFunctionReturn(PETSC_SUCCESS);
3029: }
3031: /*@C
3032: PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element
3034: Not Collective
3036: Input Parameter:
3037: + n - the number of GLL nodes
3038: . nodes - the GLL nodes
3039: . weights - the GLL weights
3041: Output Parameters:
3042: . AA - the stiffness element
3043: - AAT - the transpose of AA (pass in `NULL` if you do not need this array)
3045: Level: beginner
3047: Notes:
3048: Destroy this with `PetscGaussLobattoLegendreElementGradientDestroy()`
3050: You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented
3052: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`, `PetscGaussLobattoLegendreElementGradientDestroy()`
3053: @*/
3054: PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA, PetscReal ***AAT)
3055: {
3056: PetscReal **A, **AT = NULL;
3057: const PetscReal *gllnodes = nodes;
3058: const PetscInt p = n - 1;
3059: PetscReal Li, Lj, d0;
3060: PetscInt i, j;
3062: PetscFunctionBegin;
3063: PetscCall(PetscMalloc1(n, &A));
3064: PetscCall(PetscMalloc1(n * n, &A[0]));
3065: for (i = 1; i < n; i++) A[i] = A[i - 1] + n;
3067: if (AAT) {
3068: PetscCall(PetscMalloc1(n, &AT));
3069: PetscCall(PetscMalloc1(n * n, &AT[0]));
3070: for (i = 1; i < n; i++) AT[i] = AT[i - 1] + n;
3071: }
3073: if (n == 1) A[0][0] = 0.;
3074: d0 = (PetscReal)p * ((PetscReal)p + 1.) / 4.;
3075: for (i = 0; i < n; i++) {
3076: for (j = 0; j < n; j++) {
3077: A[i][j] = 0.;
3078: PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li));
3079: PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj));
3080: if (i != j) A[i][j] = Li / (Lj * (gllnodes[i] - gllnodes[j]));
3081: if ((j == i) && (i == 0)) A[i][j] = -d0;
3082: if (j == i && i == p) A[i][j] = d0;
3083: if (AT) AT[j][i] = A[i][j];
3084: }
3085: }
3086: if (AAT) *AAT = AT;
3087: *AA = A;
3088: PetscFunctionReturn(PETSC_SUCCESS);
3089: }
3091: /*@C
3092: PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with `PetscGaussLobattoLegendreElementGradientCreate()`
3094: Not Collective
3096: Input Parameters:
3097: + n - the number of GLL nodes
3098: . nodes - the GLL nodes
3099: . weights - the GLL weights
3100: . AA - the stiffness element
3101: - AAT - the transpose of the element
3103: Level: beginner
3105: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionCreate()`
3106: @*/
3107: PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA, PetscReal ***AAT)
3108: {
3109: PetscFunctionBegin;
3110: PetscCall(PetscFree((*AA)[0]));
3111: PetscCall(PetscFree(*AA));
3112: *AA = NULL;
3113: if (AAT) {
3114: PetscCall(PetscFree((*AAT)[0]));
3115: PetscCall(PetscFree(*AAT));
3116: *AAT = NULL;
3117: }
3118: PetscFunctionReturn(PETSC_SUCCESS);
3119: }
3121: /*@C
3122: PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element
3124: Not Collective
3126: Input Parameters:
3127: + n - the number of GLL nodes
3128: . nodes - the GLL nodes
3129: - weights - the GLL weightss
3131: Output Parameter:
3132: . AA - the stiffness element
3134: Level: beginner
3136: Notes:
3137: Destroy this with `PetscGaussLobattoLegendreElementAdvectionDestroy()`
3139: This is the same as the Gradient operator multiplied by the diagonal mass matrix
3141: You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented
3143: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionDestroy()`
3144: @*/
3145: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3146: {
3147: PetscReal **D;
3148: const PetscReal *gllweights = weights;
3149: const PetscInt glln = n;
3150: PetscInt i, j;
3152: PetscFunctionBegin;
3153: PetscCall(PetscGaussLobattoLegendreElementGradientCreate(n, nodes, weights, &D, NULL));
3154: for (i = 0; i < glln; i++) {
3155: for (j = 0; j < glln; j++) D[i][j] = gllweights[i] * D[i][j];
3156: }
3157: *AA = D;
3158: PetscFunctionReturn(PETSC_SUCCESS);
3159: }
3161: /*@C
3162: PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element created with `PetscGaussLobattoLegendreElementAdvectionCreate()`
3164: Not Collective
3166: Input Parameters:
3167: + n - the number of GLL nodes
3168: . nodes - the GLL nodes
3169: . weights - the GLL weights
3170: - A - advection
3172: Level: beginner
3174: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementAdvectionCreate()`
3175: @*/
3176: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3177: {
3178: PetscFunctionBegin;
3179: PetscCall(PetscFree((*AA)[0]));
3180: PetscCall(PetscFree(*AA));
3181: *AA = NULL;
3182: PetscFunctionReturn(PETSC_SUCCESS);
3183: }
3185: PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3186: {
3187: PetscReal **A;
3188: const PetscReal *gllweights = weights;
3189: const PetscInt glln = n;
3190: PetscInt i, j;
3192: PetscFunctionBegin;
3193: PetscCall(PetscMalloc1(glln, &A));
3194: PetscCall(PetscMalloc1(glln * glln, &A[0]));
3195: for (i = 1; i < glln; i++) A[i] = A[i - 1] + glln;
3196: if (glln == 1) A[0][0] = 0.;
3197: for (i = 0; i < glln; i++) {
3198: for (j = 0; j < glln; j++) {
3199: A[i][j] = 0.;
3200: if (j == i) A[i][j] = gllweights[i];
3201: }
3202: }
3203: *AA = A;
3204: PetscFunctionReturn(PETSC_SUCCESS);
3205: }
3207: PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3208: {
3209: PetscFunctionBegin;
3210: PetscCall(PetscFree((*AA)[0]));
3211: PetscCall(PetscFree(*AA));
3212: *AA = NULL;
3213: PetscFunctionReturn(PETSC_SUCCESS);
3214: }
3216: /*@
3217: PetscDTIndexToBary - convert an index into a barycentric coordinate.
3219: Input Parameters:
3220: + len - the desired length of the barycentric tuple (usually 1 more than the dimension it represents, so a barycentric coordinate in a triangle has length 3)
3221: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
3222: - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum)
3224: Output Parameter:
3225: . coord - will be filled with the barycentric coordinate
3227: Level: beginner
3229: Note:
3230: The indices map to barycentric coordinates in lexicographic order, where the first index is the
3231: least significant and the last index is the most significant.
3233: .seealso: `PetscDTBaryToIndex()`
3234: @*/
3235: PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[])
3236: {
3237: PetscInt c, d, s, total, subtotal, nexttotal;
3239: PetscFunctionBeginHot;
3240: PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
3241: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative");
3242: if (!len) {
3243: if (!sum && !index) PetscFunctionReturn(PETSC_SUCCESS);
3244: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
3245: }
3246: for (c = 1, total = 1; c <= len; c++) {
3247: /* total is the number of ways to have a tuple of length c with sum */
3248: if (index < total) break;
3249: total = (total * (sum + c)) / c;
3250: }
3251: PetscCheck(c <= len, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range");
3252: for (d = c; d < len; d++) coord[d] = 0;
3253: for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) {
3254: /* subtotal is the number of ways to have a tuple of length c with sum s */
3255: /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */
3256: if ((index + subtotal) >= total) {
3257: coord[--c] = sum - s;
3258: index -= (total - subtotal);
3259: sum = s;
3260: total = nexttotal;
3261: subtotal = 1;
3262: nexttotal = 1;
3263: s = 0;
3264: } else {
3265: subtotal = (subtotal * (c + s)) / (s + 1);
3266: nexttotal = (nexttotal * (c - 1 + s)) / (s + 1);
3267: s++;
3268: }
3269: }
3270: PetscFunctionReturn(PETSC_SUCCESS);
3271: }
3273: /*@
3274: PetscDTBaryToIndex - convert a barycentric coordinate to an index
3276: Input Parameters:
3277: + len - the desired length of the barycentric tuple (usually 1 more than the dimension it represents, so a barycentric coordinate in a triangle has length 3)
3278: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
3279: - coord - a barycentric coordinate with the given length and sum
3281: Output Parameter:
3282: . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum)
3284: Level: beginner
3286: Note:
3287: The indices map to barycentric coordinates in lexicographic order, where the first index is the
3288: least significant and the last index is the most significant.
3290: .seealso: `PetscDTIndexToBary`
3291: @*/
3292: PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index)
3293: {
3294: PetscInt c;
3295: PetscInt i;
3296: PetscInt total;
3298: PetscFunctionBeginHot;
3299: PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
3300: if (!len) {
3301: if (!sum) {
3302: *index = 0;
3303: PetscFunctionReturn(PETSC_SUCCESS);
3304: }
3305: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
3306: }
3307: for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c;
3308: i = total - 1;
3309: c = len - 1;
3310: sum -= coord[c];
3311: while (sum > 0) {
3312: PetscInt subtotal;
3313: PetscInt s;
3315: for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s;
3316: i -= subtotal;
3317: sum -= coord[--c];
3318: }
3319: *index = i;
3320: PetscFunctionReturn(PETSC_SUCCESS);
3321: }
3323: /*@
3324: PetscQuadratureComputePermutations - Compute permutations of quadrature points corresponding to domain orientations
3326: Input Parameter:
3327: . quad - The `PetscQuadrature`
3329: Output Parameters:
3330: + Np - The number of domain orientations
3331: - perm - An array of `IS` permutations, one for ech orientation,
3333: Level: developer
3335: .seealso: `PetscQuadratureSetCellType()`, `PetscQuadrature`
3336: @*/
3337: PetscErrorCode PetscQuadratureComputePermutations(PetscQuadrature quad, PetscInt *Np, IS *perm[])
3338: {
3339: DMPolytopeType ct;
3340: const PetscReal *xq, *wq;
3341: PetscInt dim, qdim, d, Na, o, Nq, q, qp;
3343: PetscFunctionBegin;
3344: PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &xq, &wq));
3345: PetscCall(PetscQuadratureGetCellType(quad, &ct));
3346: dim = DMPolytopeTypeGetDim(ct);
3347: Na = DMPolytopeTypeGetNumArrangments(ct);
3348: PetscCall(PetscMalloc1(Na, perm));
3349: if (Np) *Np = Na;
3350: Na /= 2;
3351: for (o = -Na; o < Na; ++o) {
3352: DM refdm;
3353: PetscInt *idx;
3354: PetscReal xi0[3] = {-1., -1., -1.}, v0[3], J[9], detJ, txq[3];
3355: PetscBool flg;
3357: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
3358: PetscCall(DMPlexOrientPoint(refdm, 0, o));
3359: PetscCall(DMPlexComputeCellGeometryFEM(refdm, 0, NULL, v0, J, NULL, &detJ));
3360: PetscCall(PetscMalloc1(Nq, &idx));
3361: for (q = 0; q < Nq; ++q) {
3362: CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q * dim], txq);
3363: for (qp = 0; qp < Nq; ++qp) {
3364: PetscReal diff = 0.;
3366: for (d = 0; d < dim; ++d) diff += PetscAbsReal(txq[d] - xq[qp * dim + d]);
3367: if (diff < PETSC_SMALL) break;
3368: }
3369: PetscCheck(qp < Nq, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Transformed quad point %" PetscInt_FMT " does not match another quad point", q);
3370: idx[q] = qp;
3371: }
3372: PetscCall(DMDestroy(&refdm));
3373: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nq, idx, PETSC_OWN_POINTER, &(*perm)[o + Na]));
3374: PetscCall(ISGetInfo((*perm)[o + Na], IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg));
3375: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ordering for orientation %" PetscInt_FMT " was not a permutation", o);
3376: PetscCall(ISSetPermutation((*perm)[o + Na]));
3377: }
3378: if (!Na) (*perm)[0] = NULL;
3379: PetscFunctionReturn(PETSC_SUCCESS);
3380: }
3382: /*@
3383: PetscDTCreateDefaultQuadrature - Create default quadrature for a given cell
3385: Not collective
3387: Input Parameters:
3388: + ct - The integration domain
3389: - qorder - The desired quadrature order
3391: Output Parameters:
3392: + q - The cell quadrature
3393: - fq - The face quadrature
3395: Level: developer
3397: .seealso: `PetscFECreateDefault()`, `PetscDTGaussTensorQuadrature()`, `PetscDTSimplexQuadrature()`, `PetscDTTensorQuadratureCreate()`
3398: @*/
3399: PetscErrorCode PetscDTCreateDefaultQuadrature(DMPolytopeType ct, PetscInt qorder, PetscQuadrature *q, PetscQuadrature *fq)
3400: {
3401: const PetscInt quadPointsPerEdge = PetscMax(qorder + 1, 1);
3402: const PetscInt dim = DMPolytopeTypeGetDim(ct);
3404: PetscFunctionBegin;
3405: switch (ct) {
3406: case DM_POLYTOPE_SEGMENT:
3407: case DM_POLYTOPE_POINT_PRISM_TENSOR:
3408: case DM_POLYTOPE_QUADRILATERAL:
3409: case DM_POLYTOPE_SEG_PRISM_TENSOR:
3410: case DM_POLYTOPE_HEXAHEDRON:
3411: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
3412: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q));
3413: PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq));
3414: break;
3415: case DM_POLYTOPE_TRIANGLE:
3416: case DM_POLYTOPE_TETRAHEDRON:
3417: PetscCall(PetscDTSimplexQuadrature(dim, 2 * qorder, PETSCDTSIMPLEXQUAD_DEFAULT, q));
3418: PetscCall(PetscDTSimplexQuadrature(dim - 1, 2 * qorder, PETSCDTSIMPLEXQUAD_DEFAULT, fq));
3419: break;
3420: case DM_POLYTOPE_TRI_PRISM:
3421: case DM_POLYTOPE_TRI_PRISM_TENSOR: {
3422: PetscQuadrature q1, q2;
3424: // TODO: this should be able to use symmetric rules, but doing so causes tests to fail
3425: PetscCall(PetscDTSimplexQuadrature(2, 2 * qorder, PETSCDTSIMPLEXQUAD_CONIC, &q1));
3426: PetscCall(PetscDTGaussTensorQuadrature(1, 1, quadPointsPerEdge, -1.0, 1.0, &q2));
3427: PetscCall(PetscDTTensorQuadratureCreate(q1, q2, q));
3428: PetscCall(PetscQuadratureDestroy(&q2));
3429: *fq = q1;
3430: /* TODO Need separate quadratures for each face */
3431: } break;
3432: default:
3433: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No quadrature for celltype %s", DMPolytopeTypes[PetscMin(ct, DM_POLYTOPE_UNKNOWN)]);
3434: }
3435: PetscFunctionReturn(PETSC_SUCCESS);
3436: }