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 <petscdmplex.h>
9: #include <petscdmshell.h>
11: #if PetscDefined(HAVE_MPFR)
12: #include <mpfr.h>
13: #endif
15: const char *const PetscDTNodeTypes_shifted[] = {"default", "gaussjacobi", "equispaced", "tanhsinh", "PetscDTNodeType", "PETSCDTNODES_", NULL};
16: const char *const *const PetscDTNodeTypes = PetscDTNodeTypes_shifted + 1;
18: const char *const PetscDTSimplexQuadratureTypes_shifted[] = {"default", "conic", "minsym", "diagsym", "PetscDTSimplexQuadratureType", "PETSCDTSIMPLEXQUAD_", NULL};
19: const char *const *const PetscDTSimplexQuadratureTypes = PetscDTSimplexQuadratureTypes_shifted + 1;
21: static PetscBool GolubWelschCite = PETSC_FALSE;
22: const char GolubWelschCitation[] = "@article{GolubWelsch1969,\n"
23: " author = {Golub and Welsch},\n"
24: " title = {Calculation of Quadrature Rules},\n"
25: " journal = {Math. Comp.},\n"
26: " volume = {23},\n"
27: " number = {106},\n"
28: " pages = {221--230},\n"
29: " year = {1969}\n}\n";
31: /* Numerical tests in src/dm/dt/tests/ex1.c show that when computing the nodes and weights of Gauss-Jacobi
32: quadrature rules:
34: - in double precision, Newton's method and Golub & Welsch both work for moderate degrees (< 100),
35: - in single precision, Newton's method starts producing incorrect roots around n = 15, but
36: the weights from Golub & Welsch become a problem before then: they produces errors
37: in computing the Jacobi-polynomial Gram matrix around n = 6.
39: So we default to Newton's method (required fewer dependencies) */
40: PetscBool PetscDTGaussQuadratureNewton_Internal = PETSC_TRUE;
42: PetscClassId PETSCQUADRATURE_CLASSID = 0;
44: /*@
45: PetscQuadratureCreate - Create a `PetscQuadrature` object
47: Collective
49: Input Parameter:
50: . comm - The communicator for the `PetscQuadrature` object
52: Output Parameter:
53: . q - The `PetscQuadrature` object
55: Level: beginner
57: .seealso: `PetscQuadrature`, `Petscquadraturedestroy()`, `PetscQuadratureGetData()`
58: @*/
59: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
60: {
61: PetscFunctionBegin;
62: PetscAssertPointer(q, 2);
63: PetscCall(DMInitializePackage());
64: PetscCall(PetscHeaderCreate(*q, PETSCQUADRATURE_CLASSID, "PetscQuadrature", "Quadrature", "DT", comm, PetscQuadratureDestroy, PetscQuadratureView));
65: (*q)->ct = DM_POLYTOPE_UNKNOWN;
66: (*q)->dim = -1;
67: (*q)->Nc = 1;
68: (*q)->order = -1;
69: (*q)->numPoints = 0;
70: (*q)->points = NULL;
71: (*q)->weights = NULL;
72: PetscFunctionReturn(PETSC_SUCCESS);
73: }
75: /*@
76: PetscQuadratureDuplicate - Create a deep copy of the `PetscQuadrature` object
78: Collective
80: Input Parameter:
81: . q - The `PetscQuadrature` object
83: Output Parameter:
84: . r - The new `PetscQuadrature` object
86: Level: beginner
88: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureDestroy()`, `PetscQuadratureGetData()`
89: @*/
90: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
91: {
92: DMPolytopeType ct;
93: PetscInt order, dim, Nc, Nq;
94: const PetscReal *points, *weights;
95: PetscReal *p, *w;
97: PetscFunctionBegin;
98: PetscAssertPointer(q, 1);
99: PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)q), r));
100: PetscCall(PetscQuadratureGetCellType(q, &ct));
101: PetscCall(PetscQuadratureSetCellType(*r, ct));
102: PetscCall(PetscQuadratureGetOrder(q, &order));
103: PetscCall(PetscQuadratureSetOrder(*r, order));
104: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights));
105: PetscCall(PetscMalloc1(Nq * dim, &p));
106: PetscCall(PetscMalloc1(Nq * Nc, &w));
107: PetscCall(PetscArraycpy(p, points, Nq * dim));
108: PetscCall(PetscArraycpy(w, weights, Nc * Nq));
109: PetscCall(PetscQuadratureSetData(*r, dim, Nc, Nq, p, w));
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: /*@
114: PetscQuadratureDestroy - Destroys a `PetscQuadrature` object
116: Collective
118: Input Parameter:
119: . q - The `PetscQuadrature` object
121: Level: beginner
123: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
124: @*/
125: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
126: {
127: PetscFunctionBegin;
128: if (!*q) PetscFunctionReturn(PETSC_SUCCESS);
130: if (--((PetscObject)*q)->refct > 0) {
131: *q = NULL;
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
134: PetscCall(PetscFree((*q)->points));
135: PetscCall(PetscFree((*q)->weights));
136: PetscCall(PetscHeaderDestroy(q));
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: /*@
141: PetscQuadratureGetCellType - Return the cell type of the integration domain
143: Not Collective
145: Input Parameter:
146: . q - The `PetscQuadrature` object
148: Output Parameter:
149: . ct - The cell type of the integration domain
151: Level: intermediate
153: .seealso: `PetscQuadrature`, `PetscQuadratureSetCellType()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
154: @*/
155: PetscErrorCode PetscQuadratureGetCellType(PetscQuadrature q, DMPolytopeType *ct)
156: {
157: PetscFunctionBegin;
159: PetscAssertPointer(ct, 2);
160: *ct = q->ct;
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /*@
165: PetscQuadratureSetCellType - Set the cell type of the integration domain
167: Not Collective
169: Input Parameters:
170: + q - The `PetscQuadrature` object
171: - ct - The cell type of the integration domain
173: Level: intermediate
175: .seealso: `PetscQuadrature`, `PetscQuadratureGetCellType()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
176: @*/
177: PetscErrorCode PetscQuadratureSetCellType(PetscQuadrature q, DMPolytopeType ct)
178: {
179: PetscFunctionBegin;
181: q->ct = ct;
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: /*@
186: PetscQuadratureGetOrder - Return the order of the method in the `PetscQuadrature`
188: Not Collective
190: Input Parameter:
191: . q - The `PetscQuadrature` object
193: Output Parameter:
194: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
196: Level: intermediate
198: .seealso: `PetscQuadrature`, `PetscQuadratureSetOrder()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
199: @*/
200: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
201: {
202: PetscFunctionBegin;
204: PetscAssertPointer(order, 2);
205: *order = q->order;
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: /*@
210: PetscQuadratureSetOrder - Set the order of the method in the `PetscQuadrature`
212: Not Collective
214: Input Parameters:
215: + q - The `PetscQuadrature` object
216: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
218: Level: intermediate
220: .seealso: `PetscQuadrature`, `PetscQuadratureGetOrder()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
221: @*/
222: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
223: {
224: PetscFunctionBegin;
226: q->order = order;
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: /*@
231: PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated
233: Not Collective
235: Input Parameter:
236: . q - The `PetscQuadrature` object
238: Output Parameter:
239: . Nc - The number of components
241: Level: intermediate
243: Note:
244: We are performing an integral $\int f(x) w(x) dx$, where both $f$ and $w$ (the weight) have `Nc` components.
246: .seealso: `PetscQuadrature`, `PetscQuadratureSetNumComponents()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
247: @*/
248: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
249: {
250: PetscFunctionBegin;
252: PetscAssertPointer(Nc, 2);
253: *Nc = q->Nc;
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: /*@
258: PetscQuadratureSetNumComponents - Sets the number of components for functions to be integrated
260: Not Collective
262: Input Parameters:
263: + q - The `PetscQuadrature` object
264: - Nc - The number of components
266: Level: intermediate
268: Note:
269: We are performing an integral $\int f(x) w(x) dx$, where both $f$ and $w$ (the weight) have `Nc` components.
271: .seealso: `PetscQuadrature`, `PetscQuadratureGetNumComponents()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
272: @*/
273: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
274: {
275: PetscFunctionBegin;
277: q->Nc = Nc;
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: /*@
282: PetscQuadratureGetData - Returns the data defining the `PetscQuadrature`
284: Not Collective
286: Input Parameter:
287: . q - The `PetscQuadrature` object
289: Output Parameters:
290: + dim - The spatial dimension
291: . Nc - The number of components
292: . npoints - The number of quadrature points
293: . points - The coordinates of each quadrature point
294: - weights - The weight of each quadrature point
296: Level: intermediate
298: Note:
299: All output arguments are optional, pass `NULL` for any argument not required
301: Fortran Note:
302: Call `PetscQuadratureRestoreData()` when you are done with the data
304: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureSetData()`
305: @*/
306: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PeOp PetscInt *dim, PeOp PetscInt *Nc, PeOp PetscInt *npoints, PeOp const PetscReal *points[], PeOp const PetscReal *weights[])
307: {
308: PetscFunctionBegin;
310: if (dim) {
311: PetscAssertPointer(dim, 2);
312: *dim = q->dim;
313: }
314: if (Nc) {
315: PetscAssertPointer(Nc, 3);
316: *Nc = q->Nc;
317: }
318: if (npoints) {
319: PetscAssertPointer(npoints, 4);
320: *npoints = q->numPoints;
321: }
322: if (points) {
323: PetscAssertPointer(points, 5);
324: *points = q->points;
325: }
326: if (weights) {
327: PetscAssertPointer(weights, 6);
328: *weights = q->weights;
329: }
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: /*@
334: PetscQuadratureEqual - determine whether two quadratures are equivalent
336: Input Parameters:
337: + A - A `PetscQuadrature` object
338: - B - Another `PetscQuadrature` object
340: Output Parameter:
341: . equal - `PETSC_TRUE` if the quadratures are the same
343: Level: intermediate
345: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`
346: @*/
347: PetscErrorCode PetscQuadratureEqual(PetscQuadrature A, PetscQuadrature B, PetscBool *equal)
348: {
349: PetscFunctionBegin;
352: PetscAssertPointer(equal, 3);
353: *equal = PETSC_FALSE;
354: if (A->ct != B->ct || A->dim != B->dim || A->Nc != B->Nc || A->order != B->order || A->numPoints != B->numPoints) PetscFunctionReturn(PETSC_SUCCESS);
355: for (PetscInt i = 0; i < A->numPoints * A->dim; i++) {
356: if (PetscAbsReal(A->points[i] - B->points[i]) > PETSC_SMALL) PetscFunctionReturn(PETSC_SUCCESS);
357: }
358: if (!A->weights && !B->weights) {
359: *equal = PETSC_TRUE;
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
362: if (A->weights && B->weights) {
363: for (PetscInt i = 0; i < A->numPoints; i++) {
364: if (PetscAbsReal(A->weights[i] - B->weights[i]) > PETSC_SMALL) PetscFunctionReturn(PETSC_SUCCESS);
365: }
366: *equal = PETSC_TRUE;
367: }
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[])
372: {
373: PetscScalar *Js, *Jinvs;
374: PetscInt i, j, k;
375: PetscBLASInt bm, bn;
377: PetscFunctionBegin;
378: if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
379: PetscCall(PetscBLASIntCast(m, &bm));
380: PetscCall(PetscBLASIntCast(n, &bn));
381: #if PetscDefined(USE_COMPLEX)
382: PetscCall(PetscMalloc2(m * n, &Js, m * n, &Jinvs));
383: for (i = 0; i < m * n; i++) Js[i] = J[i];
384: #else
385: Js = (PetscReal *)J;
386: Jinvs = Jinv;
387: #endif
388: if (m == n) {
389: PetscBLASInt *pivots;
390: PetscScalar *W;
392: PetscCall(PetscMalloc2(m, &pivots, m, &W));
394: PetscCall(PetscArraycpy(Jinvs, Js, m * m));
395: PetscCallLAPACKInfo("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info));
396: PetscCallLAPACKInfo("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &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: PetscCallLAPACKInfo("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info));
415: PetscCallLAPACKInfo("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info));
416: for (i = 0; i < n; i++) {
417: for (j = 0; j < m; j++) {
418: PetscScalar val = 0.;
420: for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j];
421: Jinvs[i * m + j] = val;
422: }
423: }
424: PetscCall(PetscFree2(pivots, W));
425: PetscCall(PetscFree(JJT));
426: } else {
427: PetscScalar *JTJ;
428: PetscBLASInt *pivots;
429: PetscScalar *W;
431: PetscCall(PetscMalloc1(n * n, &JTJ));
432: PetscCall(PetscMalloc2(n, &pivots, n, &W));
433: for (i = 0; i < n; i++) {
434: for (j = 0; j < n; j++) {
435: PetscScalar val = 0.;
437: for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j];
438: JTJ[i * n + j] = val;
439: }
440: }
442: PetscCallLAPACKInfo("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bn, pivots, &info));
443: PetscCallLAPACKInfo("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info));
444: for (i = 0; i < n; i++) {
445: for (j = 0; j < m; j++) {
446: PetscScalar val = 0.;
448: for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k];
449: Jinvs[i * m + j] = val;
450: }
451: }
452: PetscCall(PetscFree2(pivots, W));
453: PetscCall(PetscFree(JTJ));
454: }
455: #if PetscDefined(USE_COMPLEX)
456: for (i = 0; i < m * n; i++) Jinv[i] = PetscRealPart(Jinvs[i]);
457: PetscCall(PetscFree2(Js, Jinvs));
458: #endif
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: /*@
463: PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation.
465: Collective
467: Input Parameters:
468: + q - the quadrature functional
469: . imageDim - the dimension of the image of the transformation
470: . origin - a point in the original space
471: . originImage - the image of the origin under the transformation
472: . J - the Jacobian of the image: an [imageDim x dim] matrix in row major order
473: - formDegree - transform the quadrature weights as k-forms of this form degree (if the number of components is a multiple of (dim choose `formDegree`),
474: it is assumed that they represent multiple k-forms) [see `PetscDTAltVPullback()` for interpretation of `formDegree`]
476: Output Parameter:
477: . 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
478: been pulled-back by the pseudoinverse of `J` to the k-form weights in the image space.
480: Level: intermediate
482: Note:
483: The new quadrature rule will have a different number of components if spaces have different dimensions.
484: 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.
486: .seealso: `PetscQuadrature`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
487: @*/
488: PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq)
489: {
490: PetscInt dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c;
491: const PetscReal *points;
492: const PetscReal *weights;
493: PetscReal *imagePoints, *imageWeights;
494: PetscReal *Jinv;
495: PetscReal *Jinvstar;
497: PetscFunctionBegin;
499: PetscCheck(imageDim >= PetscAbsInt(formDegree), PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Cannot represent a %" PetscInt_FMT "-form in %" PetscInt_FMT " dimensions", PetscAbsInt(formDegree), imageDim);
500: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights));
501: PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize));
502: 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);
503: Ncopies = Nc / formSize;
504: PetscCall(PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize));
505: imageNc = Ncopies * imageFormSize;
506: PetscCall(PetscMalloc1(Npoints * imageDim, &imagePoints));
507: PetscCall(PetscMalloc1(Npoints * imageNc, &imageWeights));
508: PetscCall(PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar));
509: PetscCall(PetscDTJacobianInverse_Internal(imageDim, dim, J, Jinv));
510: PetscCall(PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar));
511: for (pt = 0; pt < Npoints; pt++) {
512: const PetscReal *point = PetscSafePointerPlusOffset(points, pt * dim);
513: PetscReal *imagePoint = &imagePoints[pt * imageDim];
515: for (i = 0; i < imageDim; i++) {
516: PetscReal val = originImage[i];
518: for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]);
519: imagePoint[i] = val;
520: }
521: for (c = 0; c < Ncopies; c++) {
522: const PetscReal *form = &weights[pt * Nc + c * formSize];
523: PetscReal *imageForm = &imageWeights[pt * imageNc + c * imageFormSize];
525: for (i = 0; i < imageFormSize; i++) {
526: PetscReal val = 0.;
528: for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j];
529: imageForm[i] = val;
530: }
531: }
532: }
533: PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq));
534: PetscCall(PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights));
535: PetscCall(PetscFree2(Jinv, Jinvstar));
536: PetscFunctionReturn(PETSC_SUCCESS);
537: }
539: /*@
540: PetscQuadratureSetData - Sets the data defining the quadrature
542: Not Collective
544: Input Parameters:
545: + q - The `PetscQuadrature` object
546: . dim - The spatial dimension
547: . Nc - The number of components
548: . npoints - The number of quadrature points
549: . points - The coordinates of each quadrature point
550: - weights - The weight of each quadrature point
552: Level: intermediate
554: Note:
555: `q` owns the references to points and weights, so they must be allocated using `PetscMalloc()` and the user should not free them.
557: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
558: @*/
559: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[]) PeNSS
560: {
561: PetscFunctionBegin;
563: if (dim >= 0) q->dim = dim;
564: if (Nc >= 0) q->Nc = Nc;
565: if (npoints >= 0) q->numPoints = npoints;
566: if (points) {
567: PetscAssertPointer(points, 5);
568: q->points = points;
569: }
570: if (weights) {
571: PetscAssertPointer(weights, 6);
572: q->weights = weights;
573: }
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
578: {
579: PetscInt q, d, c;
580: PetscViewerFormat format;
582: PetscFunctionBegin;
583: if (quad->Nc > 1)
584: 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));
585: 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));
586: PetscCall(PetscViewerGetFormat(v, &format));
587: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
588: for (q = 0; q < quad->numPoints; ++q) {
589: PetscCall(PetscViewerASCIIPrintf(v, "p%" PetscInt_FMT " (", q));
590: PetscCall(PetscViewerASCIIUseTabs(v, PETSC_FALSE));
591: for (d = 0; d < quad->dim; ++d) {
592: if (d) PetscCall(PetscViewerASCIIPrintf(v, ", "));
593: PetscCall(PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q * quad->dim + d]));
594: }
595: PetscCall(PetscViewerASCIIPrintf(v, ") "));
596: if (quad->Nc > 1) PetscCall(PetscViewerASCIIPrintf(v, "w%" PetscInt_FMT " (", q));
597: for (c = 0; c < quad->Nc; ++c) {
598: if (c) PetscCall(PetscViewerASCIIPrintf(v, ", "));
599: PetscCall(PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q * quad->Nc + c]));
600: }
601: if (quad->Nc > 1) PetscCall(PetscViewerASCIIPrintf(v, ")"));
602: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
603: PetscCall(PetscViewerASCIIUseTabs(v, PETSC_TRUE));
604: }
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: /*@
609: PetscQuadratureView - View a `PetscQuadrature` object
611: Collective
613: Input Parameters:
614: + quad - The `PetscQuadrature` object
615: - viewer - The `PetscViewer` object
617: Level: beginner
619: .seealso: `PetscQuadrature`, `PetscViewer`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
620: @*/
621: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
622: {
623: PetscBool isascii;
625: PetscFunctionBegin;
628: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)quad), &viewer));
629: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
630: PetscCall(PetscViewerASCIIPushTab(viewer));
631: if (isascii) PetscCall(PetscQuadratureView_Ascii(quad, viewer));
632: PetscCall(PetscViewerASCIIPopTab(viewer));
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: /*@
637: PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement
639: Not Collective; No Fortran Support
641: Input Parameters:
642: + q - The original `PetscQuadrature`
643: . numSubelements - The number of subelements the original element is divided into
644: . v0 - An array of the initial points for each subelement
645: - jac - An array of the Jacobian mappings from the reference to each subelement
647: Output Parameter:
648: . qref - The dimension
650: Level: intermediate
652: Note:
653: Together `v0` and `jac` define an affine mapping from the original reference element to each subelement
655: .seealso: `PetscQuadrature`, `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()`
656: @*/
657: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
658: {
659: DMPolytopeType ct;
660: const PetscReal *points, *weights;
661: PetscReal *pointsRef, *weightsRef;
662: PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
664: PetscFunctionBegin;
666: PetscAssertPointer(v0, 3);
667: PetscAssertPointer(jac, 4);
668: PetscAssertPointer(qref, 5);
669: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, qref));
670: PetscCall(PetscQuadratureGetCellType(q, &ct));
671: PetscCall(PetscQuadratureGetOrder(q, &order));
672: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
673: npointsRef = npoints * numSubelements;
674: PetscCall(PetscMalloc1(npointsRef * dim, &pointsRef));
675: PetscCall(PetscMalloc1(npointsRef * Nc, &weightsRef));
676: for (c = 0; c < numSubelements; ++c) {
677: for (p = 0; p < npoints; ++p) {
678: for (d = 0; d < dim; ++d) {
679: pointsRef[(c * npoints + p) * dim + d] = v0[c * dim + d];
680: for (e = 0; e < dim; ++e) pointsRef[(c * npoints + p) * dim + d] += jac[(c * dim + d) * dim + e] * (points[p * dim + e] + 1.0);
681: }
682: /* Could also use detJ here */
683: for (cp = 0; cp < Nc; ++cp) weightsRef[(c * npoints + p) * Nc + cp] = weights[p * Nc + cp] / numSubelements;
684: }
685: }
686: PetscCall(PetscQuadratureSetCellType(*qref, ct));
687: PetscCall(PetscQuadratureSetOrder(*qref, order));
688: PetscCall(PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef));
689: PetscFunctionReturn(PETSC_SUCCESS);
690: }
692: /* Compute the coefficients for the Jacobi polynomial recurrence,
694: J^{a,b}_n(x) = (cnm1 + cnm1x * x) * J^{a,b}_{n-1}(x) - cnm2 * J^{a,b}_{n-2}(x).
695: */
696: #define PetscDTJacobiRecurrence_Internal(n, a, b, cnm1, cnm1x, cnm2) \
697: do { \
698: PetscReal _a = (a); \
699: PetscReal _b = (b); \
700: PetscReal _n = (n); \
701: if (n == 1) { \
702: (cnm1) = (_a - _b) * 0.5; \
703: (cnm1x) = (_a + _b + 2.) * 0.5; \
704: (cnm2) = 0.; \
705: } else { \
706: PetscReal _2n = _n + _n; \
707: PetscReal _d = (_2n * (_n + _a + _b) * (_2n + _a + _b - 2)); \
708: PetscReal _n1 = (_2n + _a + _b - 1.) * (_a * _a - _b * _b); \
709: PetscReal _n1x = (_2n + _a + _b - 1.) * (_2n + _a + _b) * (_2n + _a + _b - 2); \
710: PetscReal _n2 = 2. * ((_n + _a - 1.) * (_n + _b - 1.) * (_2n + _a + _b)); \
711: (cnm1) = _n1 / _d; \
712: (cnm1x) = _n1x / _d; \
713: (cnm2) = _n2 / _d; \
714: } \
715: } while (0)
717: /*@
718: PetscDTJacobiNorm - Compute the weighted L2 norm of a Jacobi polynomial.
720: $\| P^{\alpha,\beta}_n \|_{\alpha,\beta}^2 = \int_{-1}^1 (1 + x)^{\alpha} (1 - x)^{\beta} P^{\alpha,\beta}_n (x)^2 dx.$
722: Input Parameters:
723: + alpha - the left exponent > -1
724: . beta - the right exponent > -1
725: - n - the polynomial degree
727: Output Parameter:
728: . norm - the weighted L2 norm
730: Level: beginner
732: .seealso: `PetscQuadrature`, `PetscDTJacobiEval()`
733: @*/
734: PetscErrorCode PetscDTJacobiNorm(PetscReal alpha, PetscReal beta, PetscInt n, PetscReal *norm)
735: {
736: PetscReal twoab1;
737: PetscReal gr;
739: PetscFunctionBegin;
740: PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent alpha %g <= -1. invalid", (double)alpha);
741: PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent beta %g <= -1. invalid", (double)beta);
742: PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "n %" PetscInt_FMT " < 0 invalid", n);
743: twoab1 = PetscPowReal(2., alpha + beta + 1.);
744: #if PetscDefined(HAVE_LGAMMA)
745: if (!n) {
746: gr = PetscExpReal(PetscLGamma(alpha + 1.) + PetscLGamma(beta + 1.) - PetscLGamma(alpha + beta + 2.));
747: } else {
748: gr = PetscExpReal(PetscLGamma(n + alpha + 1.) + PetscLGamma(n + beta + 1.) - (PetscLGamma(n + 1.) + PetscLGamma(n + alpha + beta + 1.))) / (n + n + alpha + beta + 1.);
749: }
750: #else
751: {
752: PetscInt alphai = (PetscInt)alpha;
753: PetscInt betai = (PetscInt)beta;
754: PetscInt i;
756: gr = n ? (1. / (n + n + alpha + beta + 1.)) : 1.;
757: if ((PetscReal)alphai == alpha) {
758: if (!n) {
759: for (i = 0; i < alphai; i++) gr *= (i + 1.) / (beta + i + 1.);
760: gr /= (alpha + beta + 1.);
761: } else {
762: for (i = 0; i < alphai; i++) gr *= (n + i + 1.) / (n + beta + i + 1.);
763: }
764: } else if ((PetscReal)betai == beta) {
765: if (!n) {
766: for (i = 0; i < betai; i++) gr *= (i + 1.) / (alpha + i + 2.);
767: gr /= (alpha + beta + 1.);
768: } else {
769: for (i = 0; i < betai; i++) gr *= (n + i + 1.) / (n + alpha + i + 1.);
770: }
771: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
772: }
773: #endif
774: *norm = PetscSqrtReal(twoab1 * gr);
775: PetscFunctionReturn(PETSC_SUCCESS);
776: }
778: static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p)
779: {
780: PetscReal ak, bk;
781: PetscReal abk1;
782: PetscInt i, l, maxdegree;
784: PetscFunctionBegin;
785: maxdegree = degrees[ndegree - 1] - k;
786: ak = a + k;
787: bk = b + k;
788: abk1 = a + b + k + 1.;
789: if (maxdegree < 0) {
790: for (i = 0; i < npoints; i++)
791: for (l = 0; l < ndegree; l++) p[i * ndegree + l] = 0.;
792: PetscFunctionReturn(PETSC_SUCCESS);
793: }
794: for (i = 0; i < npoints; i++) {
795: PetscReal pm1, pm2, x;
796: PetscReal cnm1, cnm1x, cnm2;
797: PetscInt j;
799: x = points[i];
800: pm2 = 1.;
801: PetscDTJacobiRecurrence_Internal(1, ak, bk, cnm1, cnm1x, cnm2);
802: pm1 = (cnm1 + cnm1x * x);
803: l = 0;
804: while (l < ndegree && degrees[l] - k < 0) p[l++] = 0.;
805: while (l < ndegree && degrees[l] - k == 0) {
806: p[l] = pm2;
807: for (PetscInt m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5;
808: l++;
809: }
810: while (l < ndegree && degrees[l] - k == 1) {
811: p[l] = pm1;
812: for (PetscInt m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5;
813: l++;
814: }
815: for (j = 2; j <= maxdegree; j++) {
816: PetscReal pp;
818: PetscDTJacobiRecurrence_Internal(j, ak, bk, cnm1, cnm1x, cnm2);
819: pp = (cnm1 + cnm1x * x) * pm1 - cnm2 * pm2;
820: pm2 = pm1;
821: pm1 = pp;
822: while (l < ndegree && degrees[l] - k == j) {
823: p[l] = pp;
824: for (PetscInt m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5;
825: l++;
826: }
827: }
828: p += ndegree;
829: }
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }
833: /*@
834: PetscDTJacobiEvalJet - Evaluate the jet (function and derivatives) of the Jacobi polynomials basis up to a given degree.
836: Input Parameters:
837: + alpha - the left exponent of the weight
838: . beta - the right exponetn of the weight
839: . npoints - the number of points to evaluate the polynomials at
840: . points - [npoints] array of point coordinates
841: . degree - the maximm degree polynomial space to evaluate, (degree + 1) will be evaluated total.
842: - k - the maximum derivative to evaluate in the jet, (k + 1) will be evaluated total.
844: Output Parameter:
845: . p - an array containing the evaluations of the Jacobi polynomials's jets on the points. the size is (degree + 1) x
846: (k + 1) x npoints, which also describes the order of the dimensions of this three-dimensional array: the first
847: (slowest varying) dimension is polynomial degree; the second dimension is derivative order; the third (fastest
848: varying) dimension is the index of the evaluation point.
850: Level: advanced
852: Notes:
853: The Jacobi polynomials with indices $\alpha$ and $\beta$ are orthogonal with respect to the
854: weighted inner product $\langle f, g \rangle = \int_{-1}^1 (1+x)^{\alpha} (1-x)^{\beta} f(x)
855: g(x) dx$.
857: .seealso: `PetscDTJacobiEval()`, `PetscDTPKDEvalJet()`
858: @*/
859: PetscErrorCode PetscDTJacobiEvalJet(PetscReal alpha, PetscReal beta, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
860: {
861: PetscInt i, j, l;
862: PetscInt *degrees;
863: PetscReal *psingle;
865: PetscFunctionBegin;
866: if (degree == 0) {
867: PetscInt zero = 0;
869: for (i = 0; i <= k; i++) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, 1, &zero, &p[i * npoints]));
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
872: PetscCall(PetscMalloc1(degree + 1, °rees));
873: PetscCall(PetscMalloc1((degree + 1) * npoints, &psingle));
874: for (i = 0; i <= degree; i++) degrees[i] = i;
875: for (i = 0; i <= k; i++) {
876: PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, degree + 1, degrees, psingle));
877: for (j = 0; j <= degree; j++) {
878: for (l = 0; l < npoints; l++) p[(j * (k + 1) + i) * npoints + l] = psingle[l * (degree + 1) + j];
879: }
880: }
881: PetscCall(PetscFree(psingle));
882: PetscCall(PetscFree(degrees));
883: PetscFunctionReturn(PETSC_SUCCESS);
884: }
886: /*@
887: PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$ at a set of points
888: at points
890: Not Collective
892: Input Parameters:
893: + npoints - number of spatial points to evaluate at
894: . alpha - the left exponent > -1
895: . beta - the right exponent > -1
896: . points - array of locations to evaluate at
897: . ndegree - number of basis degrees to evaluate
898: - degrees - sorted array of degrees to evaluate
900: Output Parameters:
901: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or `NULL`)
902: . D - row-oriented derivative evaluation matrix (or `NULL`)
903: - D2 - row-oriented second derivative evaluation matrix (or `NULL`)
905: Level: intermediate
907: .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()`
908: @*/
909: PetscErrorCode PetscDTJacobiEval(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PeOp PetscReal B[], PeOp PetscReal D[], PeOp PetscReal D2[])
910: {
911: PetscFunctionBegin;
912: PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1.");
913: PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1.");
914: if (!npoints || !ndegree) PetscFunctionReturn(PETSC_SUCCESS);
915: if (B) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B));
916: if (D) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D));
917: if (D2) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2));
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
921: /*@
922: PetscDTLegendreEval - evaluate Legendre polynomials at points
924: Not Collective
926: Input Parameters:
927: + npoints - number of spatial points to evaluate at
928: . points - array of locations to evaluate at
929: . ndegree - number of basis degrees to evaluate
930: - degrees - sorted array of degrees to evaluate
932: Output Parameters:
933: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension `npoints`*`ndegrees`, allocated by caller) (or `NULL`)
934: . D - row-oriented derivative evaluation matrix (or `NULL`)
935: - D2 - row-oriented second derivative evaluation matrix (or `NULL`)
937: Level: intermediate
939: .seealso: `PetscDTGaussQuadrature()`
940: @*/
941: PetscErrorCode PetscDTLegendreEval(PetscInt npoints, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PeOp PetscReal B[], PeOp PetscReal D[], PeOp PetscReal D2[])
942: {
943: PetscFunctionBegin;
944: PetscCall(PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2));
945: PetscFunctionReturn(PETSC_SUCCESS);
946: }
948: /*@
949: 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,
950: then the index of x is smaller than the index of y)
952: Input Parameters:
953: + len - the desired length of the degree tuple
954: - index - the index to convert: should be >= 0
956: Output Parameter:
957: . degtup - filled with a tuple of degrees
959: Level: beginner
961: Note:
962: For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
963: acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
964: last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).
966: .seealso: `PetscDTGradedOrderToIndex()`
967: @*/
968: PetscErrorCode PetscDTIndexToGradedOrder(PetscInt len, PetscInt index, PetscInt degtup[])
969: {
970: PetscInt i, total;
971: PetscInt sum;
973: PetscFunctionBeginHot;
974: PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
975: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative");
976: total = 1;
977: sum = 0;
978: while (index >= total) {
979: index -= total;
980: total = (total * (len + sum)) / (sum + 1);
981: sum++;
982: }
983: for (i = 0; i < len; i++) {
984: PetscInt c;
986: degtup[i] = sum;
987: for (c = 0, total = 1; c < sum; c++) {
988: /* going into the loop, total is the number of way to have a tuple of sum exactly c with length len - 1 - i */
989: if (index < total) break;
990: index -= total;
991: total = (total * (len - 1 - i + c)) / (c + 1);
992: degtup[i]--;
993: }
994: sum -= degtup[i];
995: }
996: PetscFunctionReturn(PETSC_SUCCESS);
997: }
999: /*@
1000: PetscDTGradedOrderToIndex - convert a tuple into an index in a graded order, the inverse of `PetscDTIndexToGradedOrder()`.
1002: Input Parameters:
1003: + len - the length of the degree tuple
1004: - degtup - tuple with this length
1006: Output Parameter:
1007: . index - index in graded order: >= 0
1009: Level: beginner
1011: Note:
1012: For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
1013: acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
1014: last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).
1016: .seealso: `PetscDTIndexToGradedOrder()`
1017: @*/
1018: PetscErrorCode PetscDTGradedOrderToIndex(PetscInt len, const PetscInt degtup[], PetscInt *index)
1019: {
1020: PetscInt i, idx, sum, total;
1022: PetscFunctionBeginHot;
1023: PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
1024: for (i = 0, sum = 0; i < len; i++) sum += degtup[i];
1025: idx = 0;
1026: total = 1;
1027: for (i = 0; i < sum; i++) {
1028: idx += total;
1029: total = (total * (len + i)) / (i + 1);
1030: }
1031: for (i = 0; i < len - 1; i++) {
1032: total = 1;
1033: sum -= degtup[i];
1034: for (PetscInt c = 0; c < sum; c++) {
1035: idx += total;
1036: total = (total * (len - 1 - i + c)) / (c + 1);
1037: }
1038: }
1039: *index = idx;
1040: PetscFunctionReturn(PETSC_SUCCESS);
1041: }
1043: static PetscBool PKDCite = PETSC_FALSE;
1044: const char PKDCitation[] = "@article{Kirby2010,\n"
1045: " title={Singularity-free evaluation of collapsed-coordinate orthogonal polynomials},\n"
1046: " author={Kirby, Robert C},\n"
1047: " journal={ACM Transactions on Mathematical Software (TOMS)},\n"
1048: " volume={37},\n"
1049: " number={1},\n"
1050: " pages={1--16},\n"
1051: " year={2010},\n"
1052: " publisher={ACM New York, NY, USA}\n}\n";
1054: /*@
1055: PetscDTPKDEvalJet - Evaluate the jet (function and derivatives) of the Proriol-Koornwinder-Dubiner (PKD) basis for
1056: the space of polynomials up to a given degree.
1058: Input Parameters:
1059: + dim - the number of variables in the multivariate polynomials
1060: . npoints - the number of points to evaluate the polynomials at
1061: . points - [npoints x dim] array of point coordinates
1062: . 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.
1063: - k - the maximum order partial derivative to evaluate in the jet. There are (dim + k choose dim) partial derivatives
1064: in the jet. Choosing k = 0 means to evaluate just the function and no derivatives
1066: Output Parameter:
1067: . p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is ((dim + degree)
1068: choose dim) x ((dim + k) choose dim) x npoints, which also describes the order of the dimensions of this
1069: three-dimensional array: the first (slowest varying) dimension is basis function index; the second dimension is jet
1070: index; the third (fastest varying) dimension is the index of the evaluation point.
1072: Level: advanced
1074: Notes:
1075: The PKD basis is L2-orthonormal on the biunit simplex (which is used as the reference element
1076: for finite elements in PETSc), which makes it a stable basis to use for evaluating
1077: polynomials in that domain.
1079: The ordering of the basis functions, and the ordering of the derivatives in the jet, both follow the graded
1080: ordering of `PetscDTIndexToGradedOrder()` and `PetscDTGradedOrderToIndex()`. For example, in 3D, the polynomial with
1081: 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);
1082: 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).
1084: The implementation uses Kirby's singularity-free evaluation algorithm, <https://doi.org/10.1145/1644001.1644006>.
1086: .seealso: `PetscDTGradedOrderToIndex()`, `PetscDTIndexToGradedOrder()`, `PetscDTJacobiEvalJet()`
1087: @*/
1088: PetscErrorCode PetscDTPKDEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
1089: {
1090: PetscInt degidx, kidx, d, pt;
1091: PetscInt Nk, Ndeg;
1092: PetscInt *ktup, *degtup;
1093: PetscReal *scales, initscale, scaleexp;
1095: PetscFunctionBegin;
1096: PetscCall(PetscCitationsRegister(PKDCitation, &PKDCite));
1097: PetscCall(PetscDTBinomialInt(dim + k, k, &Nk));
1098: PetscCall(PetscDTBinomialInt(degree + dim, degree, &Ndeg));
1099: PetscCall(PetscMalloc2(dim, °tup, dim, &ktup));
1100: PetscCall(PetscMalloc1(Ndeg, &scales));
1101: initscale = 1.;
1102: if (dim > 1) {
1103: PetscCall(PetscDTBinomial(dim, 2, &scaleexp));
1104: initscale = PetscPowReal(2., scaleexp * 0.5);
1105: }
1106: for (degidx = 0; degidx < Ndeg; degidx++) {
1107: PetscInt e;
1108: PetscInt m1idx = -1, m2idx = -1;
1109: PetscInt n;
1110: PetscInt degsum;
1111: PetscReal alpha;
1112: PetscReal cnm1, cnm1x, cnm2;
1113: PetscReal norm;
1115: PetscCall(PetscDTIndexToGradedOrder(dim, degidx, degtup));
1116: for (d = dim - 1; d >= 0; d--)
1117: if (degtup[d]) break;
1118: if (d < 0) { /* constant is 1 everywhere, all derivatives are zero */
1119: scales[degidx] = initscale;
1120: for (e = 0; e < dim; e++) {
1121: PetscCall(PetscDTJacobiNorm(e, 0., 0, &norm));
1122: scales[degidx] /= norm;
1123: }
1124: for (PetscInt i = 0; i < npoints; i++) p[degidx * Nk * npoints + i] = 1.;
1125: for (PetscInt i = 0; i < (Nk - 1) * npoints; i++) p[(degidx * Nk + 1) * npoints + i] = 0.;
1126: continue;
1127: }
1128: n = degtup[d];
1129: degtup[d]--;
1130: PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m1idx));
1131: if (degtup[d] > 0) {
1132: degtup[d]--;
1133: PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m2idx));
1134: degtup[d]++;
1135: }
1136: degtup[d]++;
1137: for (e = 0, degsum = 0; e < d; e++) degsum += degtup[e];
1138: alpha = 2 * degsum + d;
1139: PetscDTJacobiRecurrence_Internal(n, alpha, 0., cnm1, cnm1x, cnm2);
1141: scales[degidx] = initscale;
1142: for (e = 0, degsum = 0; e < dim; e++) {
1143: PetscReal ealpha;
1144: PetscReal enorm;
1146: ealpha = 2 * degsum + e;
1147: for (PetscInt f = 0; f < degsum; f++) scales[degidx] *= 2.;
1148: PetscCall(PetscDTJacobiNorm(ealpha, 0., degtup[e], &enorm));
1149: scales[degidx] /= enorm;
1150: degsum += degtup[e];
1151: }
1153: for (pt = 0; pt < npoints; pt++) {
1154: /* compute the multipliers */
1155: PetscReal thetanm1, thetanm1x, thetanm2;
1157: thetanm1x = dim - (d + 1) + 2. * points[pt * dim + d];
1158: for (e = d + 1; e < dim; e++) thetanm1x += points[pt * dim + e];
1159: thetanm1x *= 0.5;
1160: thetanm1 = (2. - (dim - (d + 1)));
1161: for (e = d + 1; e < dim; e++) thetanm1 -= points[pt * dim + e];
1162: thetanm1 *= 0.5;
1163: thetanm2 = thetanm1 * thetanm1;
1165: for (kidx = 0; kidx < Nk; kidx++) {
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 (PetscInt 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: p[(degidx * Nk + kidx) * npoints + pt] += mplty * 0.5 * (cnm1x - cnm1) * p[(m1idx * Nk + km1idx) * npoints + pt];
1183: if (m2idx >= 0) {
1184: p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm2 * thetanm1 * p[(m2idx * Nk + km1idx) * npoints + pt];
1185: /* second derivatives of -cnm2 * thetanm2 wrt x variable f,f2 is like - 0.5 * cnm2 */
1186: for (PetscInt f2 = f; f2 < dim; f2++) {
1187: PetscInt km2idx, mplty2 = ktup[f2];
1188: PetscInt factor;
1190: if (!mplty2) continue;
1191: ktup[f2]--;
1192: PetscCall(PetscDTGradedOrderToIndex(dim, ktup, &km2idx));
1194: factor = mplty * mplty2;
1195: if (f == f2) factor /= 2;
1196: p[(degidx * Nk + kidx) * npoints + pt] -= 0.5 * factor * cnm2 * p[(m2idx * Nk + km2idx) * npoints + pt];
1197: ktup[f2]++;
1198: }
1199: }
1200: } else {
1201: p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm1x * p[(m1idx * Nk + km1idx) * npoints + pt];
1202: }
1203: ktup[f]++;
1204: }
1205: }
1206: }
1207: }
1208: for (degidx = 0; degidx < Ndeg; degidx++) {
1209: PetscReal scale = scales[degidx];
1211: for (PetscInt i = 0; i < Nk * npoints; i++) p[degidx * Nk * npoints + i] *= scale;
1212: }
1213: PetscCall(PetscFree(scales));
1214: PetscCall(PetscFree2(degtup, ktup));
1215: PetscFunctionReturn(PETSC_SUCCESS);
1216: }
1218: /*@
1219: PetscDTPTrimmedSize - The size of the trimmed polynomial space of k-forms with a given degree and form degree,
1220: which can be evaluated in `PetscDTPTrimmedEvalJet()`.
1222: Input Parameters:
1223: + dim - the number of variables in the multivariate polynomials
1224: . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space.
1225: - formDegree - the degree of the form
1227: Output Parameter:
1228: . size - The number ((`dim` + `degree`) choose (`dim` + `formDegree`)) x ((`degree` + `formDegree` - 1) choose (`formDegree`))
1230: Level: advanced
1232: .seealso: `PetscDTPTrimmedEvalJet()`
1233: @*/
1234: PetscErrorCode PetscDTPTrimmedSize(PetscInt dim, PetscInt degree, PetscInt formDegree, PetscInt *size)
1235: {
1236: PetscInt Nrk, Nbpt; // number of trimmed polynomials
1238: PetscFunctionBegin;
1239: formDegree = PetscAbsInt(formDegree);
1240: PetscCall(PetscDTBinomialInt(degree + dim, degree + formDegree, &Nbpt));
1241: PetscCall(PetscDTBinomialInt(degree + formDegree - 1, formDegree, &Nrk));
1242: Nbpt *= Nrk;
1243: *size = Nbpt;
1244: PetscFunctionReturn(PETSC_SUCCESS);
1245: }
1247: /* there was a reference implementation based on section 4.4 of Arnold, Falk & Winther (acta numerica, 2006), but it
1248: * was inferior to this implementation */
1249: static PetscErrorCode PetscDTPTrimmedEvalJet_Internal(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[])
1250: {
1251: PetscInt formDegreeOrig = formDegree;
1252: PetscBool formNegative = (formDegreeOrig < 0) ? PETSC_TRUE : PETSC_FALSE;
1254: PetscFunctionBegin;
1255: formDegree = PetscAbsInt(formDegreeOrig);
1256: if (formDegree == 0) {
1257: PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree, jetDegree, p));
1258: PetscFunctionReturn(PETSC_SUCCESS);
1259: }
1260: if (formDegree == dim) {
1261: PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p));
1262: PetscFunctionReturn(PETSC_SUCCESS);
1263: }
1264: PetscInt Nbpt;
1265: PetscCall(PetscDTPTrimmedSize(dim, degree, formDegree, &Nbpt));
1266: PetscInt Nf;
1267: PetscCall(PetscDTBinomialInt(dim, formDegree, &Nf));
1268: PetscInt Nk;
1269: PetscCall(PetscDTBinomialInt(dim + jetDegree, dim, &Nk));
1270: PetscCall(PetscArrayzero(p, Nbpt * Nf * Nk * npoints));
1272: PetscInt Nbpm1; // number of scalar polynomials up to degree - 1;
1273: PetscCall(PetscDTBinomialInt(dim + degree - 1, dim, &Nbpm1));
1274: PetscReal *p_scalar;
1275: PetscCall(PetscMalloc1(Nbpm1 * Nk * npoints, &p_scalar));
1276: PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p_scalar));
1277: PetscInt total = 0;
1278: // First add the full polynomials up to degree - 1 into the basis: take the scalar
1279: // and copy one for each form component
1280: for (PetscInt i = 0; i < Nbpm1; i++) {
1281: const PetscReal *src = &p_scalar[i * Nk * npoints];
1282: for (PetscInt f = 0; f < Nf; f++) {
1283: PetscReal *dest = &p[(total++ * Nf + f) * Nk * npoints];
1284: PetscCall(PetscArraycpy(dest, src, Nk * npoints));
1285: }
1286: }
1287: PetscInt *form_atoms;
1288: PetscCall(PetscMalloc1(formDegree + 1, &form_atoms));
1289: // construct the interior product pattern
1290: PetscInt (*pattern)[3];
1291: PetscInt Nf1; // number of formDegree + 1 forms
1292: PetscCall(PetscDTBinomialInt(dim, formDegree + 1, &Nf1));
1293: PetscInt nnz = Nf1 * (formDegree + 1);
1294: PetscCall(PetscMalloc1(Nf1 * (formDegree + 1), &pattern));
1295: PetscCall(PetscDTAltVInteriorPattern(dim, formDegree + 1, pattern));
1296: PetscReal centroid = (1. - dim) / (dim + 1.);
1297: PetscInt *deriv;
1298: PetscCall(PetscMalloc1(dim, &deriv));
1299: for (PetscInt d = dim; d >= formDegree + 1; d--) {
1300: PetscInt Nfd1; // number of formDegree + 1 forms in dimension d that include dx_0
1301: // (equal to the number of formDegree forms in dimension d-1)
1302: PetscCall(PetscDTBinomialInt(d - 1, formDegree, &Nfd1));
1303: // The number of homogeneous (degree-1) scalar polynomials in d variables
1304: PetscInt Nh;
1305: PetscCall(PetscDTBinomialInt(d - 1 + degree - 1, d - 1, &Nh));
1306: const PetscReal *h_scalar = &p_scalar[(Nbpm1 - Nh) * Nk * npoints];
1307: for (PetscInt b = 0; b < Nh; b++) {
1308: const PetscReal *h_s = &h_scalar[b * Nk * npoints];
1309: for (PetscInt f = 0; f < Nfd1; f++) {
1310: // construct all formDegree+1 forms that start with dx_(dim - d) /\ ...
1311: form_atoms[0] = dim - d;
1312: PetscCall(PetscDTEnumSubset(d - 1, formDegree, f, &form_atoms[1]));
1313: for (PetscInt i = 0; i < formDegree; i++) form_atoms[1 + i] += form_atoms[0] + 1;
1314: PetscInt f_ind; // index of the resulting form
1315: PetscCall(PetscDTSubsetIndex(dim, formDegree + 1, form_atoms, &f_ind));
1316: PetscReal *p_f = &p[total++ * Nf * Nk * npoints];
1317: for (PetscInt nz = 0; nz < nnz; nz++) {
1318: PetscInt i = pattern[nz][0]; // formDegree component
1319: PetscInt j = pattern[nz][1]; // (formDegree + 1) component
1320: PetscInt v = pattern[nz][2]; // coordinate component
1321: PetscReal scale = v < 0 ? -1. : 1.;
1323: i = formNegative ? (Nf - 1 - i) : i;
1324: scale = (formNegative && (i & 1)) ? -scale : scale;
1325: v = v < 0 ? -(v + 1) : v;
1326: if (j != f_ind) continue;
1327: PetscReal *p_i = &p_f[i * Nk * npoints];
1328: for (PetscInt jet = 0; jet < Nk; jet++) {
1329: const PetscReal *h_jet = &h_s[jet * npoints];
1330: PetscReal *p_jet = &p_i[jet * npoints];
1332: for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * h_jet[pt] * (points[pt * dim + v] - centroid);
1333: PetscCall(PetscDTIndexToGradedOrder(dim, jet, deriv));
1334: deriv[v]++;
1335: PetscReal mult = deriv[v];
1336: PetscInt l;
1337: PetscCall(PetscDTGradedOrderToIndex(dim, deriv, &l));
1338: if (l >= Nk) continue;
1339: p_jet = &p_i[l * npoints];
1340: for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * mult * h_jet[pt];
1341: deriv[v]--;
1342: }
1343: }
1344: }
1345: }
1346: }
1347: PetscCheck(total == Nbpt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrectly counted P trimmed polynomials");
1348: PetscCall(PetscFree(deriv));
1349: PetscCall(PetscFree(pattern));
1350: PetscCall(PetscFree(form_atoms));
1351: PetscCall(PetscFree(p_scalar));
1352: PetscFunctionReturn(PETSC_SUCCESS);
1353: }
1355: /*@
1356: PetscDTPTrimmedEvalJet - Evaluate the jet (function and derivatives) of a basis of the trimmed polynomial k-forms up to
1357: a given degree.
1359: Input Parameters:
1360: + dim - the number of variables in the multivariate polynomials
1361: . npoints - the number of points to evaluate the polynomials at
1362: . points - [npoints x dim] array of point coordinates
1363: . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space to evaluate.
1364: There are ((dim + degree) choose (dim + formDegree)) x ((degree + formDegree - 1) choose (formDegree)) polynomials in this space.
1365: (You can use `PetscDTPTrimmedSize()` to compute this size.)
1366: . formDegree - the degree of the form
1367: - jetDegree - the maximum order partial derivative to evaluate in the jet. There are ((dim + jetDegree) choose dim) partial derivatives
1368: in the jet. Choosing jetDegree = 0 means to evaluate just the function and no derivatives
1370: Output Parameter:
1371: . p - an array containing the evaluations of the PKD polynomials' jets on the points.
1373: Level: advanced
1375: Notes:
1376: The size of `p` is `PetscDTPTrimmedSize()` x ((dim + formDegree) choose dim) x ((dim + k)
1377: choose dim) x npoints,which also describes the order of the dimensions of this
1378: four-dimensional array\:
1380: the first (slowest varying) dimension is basis function index;
1381: the second dimension is component of the form;
1382: the third dimension is jet index;
1383: the fourth (fastest varying) dimension is the index of the evaluation point.
1385: The ordering of the basis functions is not graded, so the basis functions are not nested by degree like `PetscDTPKDEvalJet()`.
1386: The basis functions are not an L2-orthonormal basis on any particular domain.
1388: The implementation is based on the description of the trimmed polynomials up to degree r as
1389: the direct sum of polynomials up to degree (r-1) and the Koszul differential applied to
1390: homogeneous polynomials of degree (r-1).
1392: .seealso: `PetscDTPKDEvalJet()`, `PetscDTPTrimmedSize()`
1393: @*/
1394: PetscErrorCode PetscDTPTrimmedEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[])
1395: {
1396: PetscFunctionBegin;
1397: PetscCall(PetscDTPTrimmedEvalJet_Internal(dim, npoints, points, degree, formDegree, jetDegree, p));
1398: PetscFunctionReturn(PETSC_SUCCESS);
1399: }
1401: /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V
1402: * with lds n; diag and subdiag are overwritten */
1403: static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[], PetscReal eigs[], PetscScalar V[])
1404: {
1405: char jobz = 'V'; /* eigenvalues and eigenvectors */
1406: char range = 'A'; /* all eigenvalues will be found */
1407: PetscReal VL = 0.; /* ignored because range is 'A' */
1408: PetscReal VU = 0.; /* ignored because range is 'A' */
1409: PetscBLASInt IL = 0; /* ignored because range is 'A' */
1410: PetscBLASInt IU = 0; /* ignored because range is 'A' */
1411: PetscReal abstol = 0.; /* unused */
1412: PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */
1413: PetscBLASInt *isuppz;
1414: PetscBLASInt lwork, liwork;
1415: PetscReal workquery;
1416: PetscBLASInt iworkquery;
1417: PetscBLASInt *iwork;
1418: PetscReal *work = NULL;
1420: PetscFunctionBegin;
1421: #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1422: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1423: #endif
1424: PetscCall(PetscBLASIntCast(n, &bn));
1425: PetscCall(PetscBLASIntCast(n, &ldz));
1426: #if !PetscDefined(MISSING_LAPACK_STEGR)
1427: PetscCall(PetscMalloc1(2 * n, &isuppz));
1428: lwork = -1;
1429: liwork = -1;
1430: PetscCallLAPACKInfo("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, &workquery, &lwork, &iworkquery, &liwork, &info));
1431: lwork = (PetscBLASInt)workquery;
1432: liwork = iworkquery;
1433: PetscCall(PetscMalloc2(lwork, &work, liwork, &iwork));
1434: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1435: PetscCallLAPACKInfo("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, work, &lwork, iwork, &liwork, &info));
1436: PetscCall(PetscFPTrapPop());
1437: PetscCall(PetscFree2(work, iwork));
1438: PetscCall(PetscFree(isuppz));
1439: #elif !PetscDefined(MISSING_LAPACK_STEQR)
1440: jobz = 'I'; /* Compute eigenvalues and eigenvectors of the
1441: tridiagonal matrix. Z is initialized to the identity
1442: matrix. */
1443: PetscCall(PetscMalloc1(PetscMax(1, 2 * n - 2), &work));
1444: PetscCallLAPACKInfo("LAPACKsteqr", LAPACKsteqr_("I", &bn, diag, subdiag, V, &ldz, work, &info));
1445: PetscCall(PetscFPTrapPop());
1446: PetscCall(PetscFree(work));
1447: PetscCall(PetscArraycpy(eigs, diag, n));
1448: #endif
1449: PetscFunctionReturn(PETSC_SUCCESS);
1450: }
1452: /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi
1453: * quadrature rules on the interval [-1, 1] */
1454: static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw)
1455: {
1456: PetscReal twoab1;
1457: PetscInt m = n - 2;
1458: PetscReal a = alpha + 1.;
1459: PetscReal b = beta + 1.;
1460: PetscReal gra, grb;
1462: PetscFunctionBegin;
1463: twoab1 = PetscPowReal(2., a + b - 1.);
1464: #if PetscDefined(HAVE_LGAMMA)
1465: grb = PetscExpReal(2. * PetscLGamma(b + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + a + 1.) - (PetscLGamma(m + b + 1) + PetscLGamma(m + a + b + 1.)));
1466: gra = PetscExpReal(2. * PetscLGamma(a + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + b + 1.) - (PetscLGamma(m + a + 1) + PetscLGamma(m + a + b + 1.)));
1467: #else
1468: {
1469: PetscReal binom1, binom2;
1470: PetscInt alphai = (PetscInt)alpha;
1471: PetscInt betai = (PetscInt)beta;
1473: PetscCheck((PetscReal)alphai == alpha && (PetscReal)betai == beta, PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
1474: PetscCall(PetscDTBinomial(m + b, b, &binom1));
1475: PetscCall(PetscDTBinomial(m + a + b, b, &binom2));
1476: grb = 1. / (binom1 * binom2);
1477: PetscCall(PetscDTBinomial(m + a, a, &binom1));
1478: PetscCall(PetscDTBinomial(m + a + b, a, &binom2));
1479: gra = 1. / (binom1 * binom2);
1480: }
1481: #endif
1482: *leftw = twoab1 * grb / b;
1483: *rightw = twoab1 * gra / a;
1484: PetscFunctionReturn(PETSC_SUCCESS);
1485: }
1487: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
1488: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
1489: static inline PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
1490: {
1491: PetscReal pn1, pn2;
1492: PetscReal cnm1, cnm1x, cnm2;
1494: PetscFunctionBegin;
1495: if (!n) {
1496: *P = 1.0;
1497: PetscFunctionReturn(PETSC_SUCCESS);
1498: }
1499: PetscDTJacobiRecurrence_Internal(1, a, b, cnm1, cnm1x, cnm2);
1500: pn2 = 1.;
1501: pn1 = cnm1 + cnm1x * x;
1502: if (n == 1) {
1503: *P = pn1;
1504: PetscFunctionReturn(PETSC_SUCCESS);
1505: }
1506: *P = 0.0;
1507: for (PetscInt k = 2; k < n + 1; ++k) {
1508: PetscDTJacobiRecurrence_Internal(k, a, b, cnm1, cnm1x, cnm2);
1510: *P = (cnm1 + cnm1x * x) * pn1 - cnm2 * pn2;
1511: pn2 = pn1;
1512: pn1 = *P;
1513: }
1514: PetscFunctionReturn(PETSC_SUCCESS);
1515: }
1517: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
1518: static inline PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P)
1519: {
1520: PetscReal nP;
1521: PetscInt i;
1523: PetscFunctionBegin;
1524: *P = 0.0;
1525: if (k > n) PetscFunctionReturn(PETSC_SUCCESS);
1526: PetscCall(PetscDTComputeJacobi(a + k, b + k, n - k, x, &nP));
1527: for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5;
1528: *P = nP;
1529: PetscFunctionReturn(PETSC_SUCCESS);
1530: }
1532: static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[])
1533: {
1534: PetscInt maxIter = 100;
1535: PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON));
1536: PetscReal a1, a6, gf;
1538: PetscFunctionBegin;
1539: a1 = PetscPowReal(2.0, a + b + 1);
1540: #if PetscDefined(HAVE_LGAMMA)
1541: {
1542: PetscReal a2, a3, a4, a5;
1543: a2 = PetscLGamma(a + npoints + 1);
1544: a3 = PetscLGamma(b + npoints + 1);
1545: a4 = PetscLGamma(a + b + npoints + 1);
1546: a5 = PetscLGamma(npoints + 1);
1547: gf = PetscExpReal(a2 + a3 - (a4 + a5));
1548: }
1549: #else
1550: {
1551: PetscInt ia, ib;
1553: ia = (PetscInt)a;
1554: ib = (PetscInt)b;
1555: gf = 1.;
1556: if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */
1557: for (PetscInt k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k);
1558: } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */
1559: for (PetscInt k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k);
1560: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
1561: }
1562: #endif
1564: a6 = a1 * gf;
1565: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
1566: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
1567: for (PetscInt k = 0; k < npoints; ++k) {
1568: PetscReal r = PetscCosReal(PETSC_PI * (1. - (4. * k + 3. + 2. * b) / (4. * npoints + 2. * (a + b + 1.)))), dP;
1570: if (k > 0) r = 0.5 * (r + x[k - 1]);
1571: for (PetscInt j = 0; j < maxIter; ++j) {
1572: PetscReal s = 0.0, delta, f, fp;
1573: PetscInt i;
1575: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
1576: PetscCall(PetscDTComputeJacobi(a, b, npoints, r, &f));
1577: PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp));
1578: delta = f / (fp - f * s);
1579: r = r - delta;
1580: if (PetscAbsReal(delta) < eps) break;
1581: }
1582: x[k] = r;
1583: PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP));
1584: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
1585: }
1586: PetscFunctionReturn(PETSC_SUCCESS);
1587: }
1589: /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi
1590: * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */
1591: static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s)
1592: {
1593: PetscInt i;
1595: PetscFunctionBegin;
1596: for (i = 0; i < nPoints; i++) {
1597: PetscReal A, B, C;
1599: PetscDTJacobiRecurrence_Internal(i + 1, a, b, A, B, C);
1600: d[i] = -A / B;
1601: if (i) s[i - 1] *= C / B;
1602: if (i < nPoints - 1) s[i] = 1. / B;
1603: }
1604: PetscFunctionReturn(PETSC_SUCCESS);
1605: }
1607: static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
1608: {
1609: PetscReal mu0;
1610: PetscReal ga, gb, gab;
1611: PetscInt i;
1613: PetscFunctionBegin;
1614: PetscCall(PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite));
1616: #if PetscDefined(HAVE_TGAMMA)
1617: ga = PetscTGamma(a + 1);
1618: gb = PetscTGamma(b + 1);
1619: gab = PetscTGamma(a + b + 2);
1620: #else
1621: {
1622: PetscInt ia = (PetscInt)a, ib = (PetscInt)b;
1624: PetscCheck(ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "tgamma() - math routine is unavailable.");
1625: /* All gamma(x) terms are (x-1)! terms */
1626: PetscCall(PetscDTFactorial(ia, &ga));
1627: PetscCall(PetscDTFactorial(ib, &gb));
1628: PetscCall(PetscDTFactorial(ia + ib + 1, &gab));
1629: }
1630: #endif
1631: mu0 = PetscPowReal(2., a + b + 1.) * ga * gb / gab;
1633: #if defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1634: {
1635: PetscReal *diag, *subdiag;
1636: PetscScalar *V;
1638: PetscCall(PetscMalloc2(npoints, &diag, npoints, &subdiag));
1639: PetscCall(PetscMalloc1(npoints * npoints, &V));
1640: PetscCall(PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag));
1641: for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]);
1642: PetscCall(PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V));
1643: for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0;
1644: PetscCall(PetscFree(V));
1645: PetscCall(PetscFree2(diag, subdiag));
1646: }
1647: #else
1648: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1649: #endif
1650: { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the
1651: eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that
1652: the eigenvalues are sorted */
1653: PetscBool sorted;
1655: PetscCall(PetscSortedReal(npoints, x, &sorted));
1656: if (!sorted) {
1657: PetscInt *order, i;
1658: PetscReal *tmp;
1660: PetscCall(PetscMalloc2(npoints, &order, npoints, &tmp));
1661: for (i = 0; i < npoints; i++) order[i] = i;
1662: PetscCall(PetscSortRealWithPermutation(npoints, x, order));
1663: PetscCall(PetscArraycpy(tmp, x, npoints));
1664: for (i = 0; i < npoints; i++) x[i] = tmp[order[i]];
1665: PetscCall(PetscArraycpy(tmp, w, npoints));
1666: for (i = 0; i < npoints; i++) w[i] = tmp[order[i]];
1667: PetscCall(PetscFree2(order, tmp));
1668: }
1669: }
1670: PetscFunctionReturn(PETSC_SUCCESS);
1671: }
1673: static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1674: {
1675: PetscFunctionBegin;
1676: PetscCheck(npoints >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive");
1677: /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1678: PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1.");
1679: PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1.");
1681: if (newton) PetscCall(PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w));
1682: else PetscCall(PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w));
1683: if (alpha == beta) { /* symmetrize */
1684: PetscInt i;
1685: for (i = 0; i < (npoints + 1) / 2; i++) {
1686: PetscInt j = npoints - 1 - i;
1687: PetscReal xi = x[i];
1688: PetscReal xj = x[j];
1689: PetscReal wi = w[i];
1690: PetscReal wj = w[j];
1692: x[i] = (xi - xj) / 2.;
1693: x[j] = (xj - xi) / 2.;
1694: w[i] = w[j] = (wi + wj) / 2.;
1695: }
1696: }
1697: PetscFunctionReturn(PETSC_SUCCESS);
1698: }
1700: /*@
1701: PetscDTGaussJacobiQuadrature - quadrature for the interval $[a, b]$ with the weight function
1702: $(x-a)^\alpha (x-b)^\beta$.
1704: Not Collective
1706: Input Parameters:
1707: + npoints - the number of points in the quadrature rule
1708: . a - the left endpoint of the interval
1709: . b - the right endpoint of the interval
1710: . alpha - the left exponent
1711: - beta - the right exponent
1713: Output Parameters:
1714: + x - array of length `npoints`, the locations of the quadrature points
1715: - w - array of length `npoints`, the weights of the quadrature points
1717: Level: intermediate
1719: Note:
1720: This quadrature rule is exact for polynomials up to degree 2*`npoints` - 1.
1722: .seealso: `PetscDTGaussQuadrature()`
1723: @*/
1724: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1725: {
1726: PetscInt i;
1728: PetscFunctionBegin;
1729: PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal));
1730: if (a != -1. || b != 1.) { /* shift */
1731: for (i = 0; i < npoints; i++) {
1732: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1733: w[i] *= (b - a) / 2.;
1734: }
1735: }
1736: PetscFunctionReturn(PETSC_SUCCESS);
1737: }
1739: static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1740: {
1741: PetscFunctionBegin;
1742: PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive");
1743: /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1744: PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1.");
1745: PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1.");
1747: x[0] = -1.;
1748: x[npoints - 1] = 1.;
1749: if (npoints > 2) PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints - 2, alpha + 1., beta + 1., &x[1], &w[1], newton));
1750: for (PetscInt i = 1; i < npoints - 1; i++) w[i] /= (1. - x[i] * x[i]);
1751: PetscCall(PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints - 1]));
1752: PetscFunctionReturn(PETSC_SUCCESS);
1753: }
1755: /*@
1756: PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval $[a, b]$ with the weight function
1757: $(x-a)^\alpha (x-b)^\beta$, with endpoints `a` and `b` included as quadrature points.
1759: Not Collective
1761: Input Parameters:
1762: + npoints - the number of points in the quadrature rule
1763: . a - the left endpoint of the interval
1764: . b - the right endpoint of the interval
1765: . alpha - the left exponent
1766: - beta - the right exponent
1768: Output Parameters:
1769: + x - array of length `npoints`, the locations of the quadrature points
1770: - w - array of length `npoints`, the weights of the quadrature points
1772: Level: intermediate
1774: Note:
1775: This quadrature rule is exact for polynomials up to degree 2*`npoints` - 3.
1777: .seealso: `PetscDTGaussJacobiQuadrature()`
1778: @*/
1779: PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1780: {
1781: PetscFunctionBegin;
1782: PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal));
1783: if (a != -1. || b != 1.) { /* shift */
1784: for (PetscInt i = 0; i < npoints; i++) {
1785: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1786: w[i] *= (b - a) / 2.;
1787: }
1788: }
1789: PetscFunctionReturn(PETSC_SUCCESS);
1790: }
1792: /*@
1793: PetscDTGaussQuadrature - create Gauss-Legendre quadrature
1795: Not Collective
1797: Input Parameters:
1798: + npoints - number of points
1799: . a - left end of interval (often-1)
1800: - b - right end of interval (often +1)
1802: Output Parameters:
1803: + x - quadrature points
1804: - w - quadrature weights
1806: Level: intermediate
1808: Note:
1809: See {cite}`golub1969calculation`
1811: .seealso: `PetscDTLegendreEval()`, `PetscDTGaussJacobiQuadrature()`
1812: @*/
1813: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[])
1814: {
1815: PetscFunctionBegin;
1816: PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal));
1817: if (a != -1. || b != 1.) { /* shift */
1818: for (PetscInt i = 0; i < npoints; i++) {
1819: x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1820: w[i] *= (b - a) / 2.;
1821: }
1822: }
1823: PetscFunctionReturn(PETSC_SUCCESS);
1824: }
1826: /*@
1827: PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre
1828: nodes of a given size on the domain $[-1,1]$
1830: Not Collective
1832: Input Parameters:
1833: + npoints - number of grid nodes
1834: - type - `PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA` or `PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON`
1836: Output Parameters:
1837: + x - quadrature points, pass in an array of length `npoints`
1838: - w - quadrature weights, pass in an array of length `npoints`
1840: Level: intermediate
1842: Notes:
1843: For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
1844: close enough to the desired solution
1846: These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes
1848: 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
1850: .seealso: `PetscDTGaussQuadrature()`, `PetscGaussLobattoLegendreCreateType`
1851: @*/
1852: PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints, PetscGaussLobattoLegendreCreateType type, PetscReal x[], PetscReal w[])
1853: {
1854: PetscBool newton;
1856: PetscFunctionBegin;
1857: PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide at least 2 grid points per element");
1858: newton = (PetscBool)(type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON);
1859: PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton));
1860: PetscFunctionReturn(PETSC_SUCCESS);
1861: }
1863: /*@
1864: PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
1866: Not Collective
1868: Input Parameters:
1869: + dim - The spatial dimension
1870: . Nc - The number of components
1871: . npoints - number of points in one dimension
1872: . a - left end of interval (often-1)
1873: - b - right end of interval (often +1)
1875: Output Parameter:
1876: . q - A `PetscQuadrature` object
1878: Level: intermediate
1880: .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()`
1881: @*/
1882: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1883: {
1884: DMPolytopeType ct;
1885: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints * PetscSqr(npoints) : PetscSqr(npoints) : npoints;
1886: PetscReal *x, *w, *xw, *ww;
1888: PetscFunctionBegin;
1889: PetscCall(PetscMalloc1(totpoints * dim, &x));
1890: PetscCall(PetscMalloc1(totpoints * Nc, &w));
1891: /* Set up the Golub-Welsch system */
1892: switch (dim) {
1893: case 0:
1894: ct = DM_POLYTOPE_POINT;
1895: PetscCall(PetscFree(x));
1896: PetscCall(PetscFree(w));
1897: PetscCall(PetscMalloc1(1, &x));
1898: PetscCall(PetscMalloc1(Nc, &w));
1899: totpoints = 1;
1900: x[0] = 0.0;
1901: for (PetscInt c = 0; c < Nc; ++c) w[c] = 1.0;
1902: break;
1903: case 1:
1904: ct = DM_POLYTOPE_SEGMENT;
1905: PetscCall(PetscMalloc1(npoints, &ww));
1906: PetscCall(PetscDTGaussQuadrature(npoints, a, b, x, ww));
1907: for (PetscInt i = 0; i < npoints; ++i)
1908: for (PetscInt c = 0; c < Nc; ++c) w[i * Nc + c] = ww[i];
1909: PetscCall(PetscFree(ww));
1910: break;
1911: case 2:
1912: ct = DM_POLYTOPE_QUADRILATERAL;
1913: PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww));
1914: PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww));
1915: for (PetscInt i = 0; i < npoints; ++i) {
1916: for (PetscInt j = 0; j < npoints; ++j) {
1917: x[(i * npoints + j) * dim + 0] = xw[i];
1918: x[(i * npoints + j) * dim + 1] = xw[j];
1919: for (PetscInt c = 0; c < Nc; ++c) w[(i * npoints + j) * Nc + c] = ww[i] * ww[j];
1920: }
1921: }
1922: PetscCall(PetscFree2(xw, ww));
1923: break;
1924: case 3:
1925: ct = DM_POLYTOPE_HEXAHEDRON;
1926: PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww));
1927: PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww));
1928: for (PetscInt i = 0; i < npoints; ++i) {
1929: for (PetscInt j = 0; j < npoints; ++j) {
1930: for (PetscInt k = 0; k < npoints; ++k) {
1931: x[((i * npoints + j) * npoints + k) * dim + 0] = xw[i];
1932: x[((i * npoints + j) * npoints + k) * dim + 1] = xw[j];
1933: x[((i * npoints + j) * npoints + k) * dim + 2] = xw[k];
1934: for (PetscInt c = 0; c < Nc; ++c) w[((i * npoints + j) * npoints + k) * Nc + c] = ww[i] * ww[j] * ww[k];
1935: }
1936: }
1937: }
1938: PetscCall(PetscFree2(xw, ww));
1939: break;
1940: default:
1941: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %" PetscInt_FMT, dim);
1942: }
1943: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
1944: PetscCall(PetscQuadratureSetCellType(*q, ct));
1945: PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1));
1946: PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w));
1947: PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "GaussTensor"));
1948: PetscFunctionReturn(PETSC_SUCCESS);
1949: }
1951: /*@
1952: PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex {cite}`karniadakis2005spectral`
1954: Not Collective
1956: Input Parameters:
1957: + dim - The simplex dimension
1958: . Nc - The number of components
1959: . npoints - The number of points in one dimension
1960: . a - left end of interval (often-1)
1961: - b - right end of interval (often +1)
1963: Output Parameter:
1964: . q - A `PetscQuadrature` object
1966: Level: intermediate
1968: Note:
1969: For `dim` == 1, this is Gauss-Legendre quadrature
1971: .seealso: `PetscDTGaussTensorQuadrature()`, `PetscDTGaussQuadrature()`
1972: @*/
1973: PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1974: {
1975: DMPolytopeType ct;
1976: PetscInt totpoints;
1977: PetscReal *p1, *w1;
1978: PetscReal *x, *w;
1980: PetscFunctionBegin;
1981: PetscCheck(!(a != -1.0) && !(b != 1.0), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
1982: switch (dim) {
1983: case 0:
1984: ct = DM_POLYTOPE_POINT;
1985: break;
1986: case 1:
1987: ct = DM_POLYTOPE_SEGMENT;
1988: break;
1989: case 2:
1990: ct = DM_POLYTOPE_TRIANGLE;
1991: break;
1992: case 3:
1993: ct = DM_POLYTOPE_TETRAHEDRON;
1994: break;
1995: default:
1996: ct = DM_POLYTOPE_UNKNOWN;
1997: }
1998: totpoints = 1;
1999: for (PetscInt i = 0; i < dim; ++i) totpoints *= npoints;
2000: PetscCall(PetscMalloc1(totpoints * dim, &x));
2001: PetscCall(PetscMalloc1(totpoints * Nc, &w));
2002: PetscCall(PetscMalloc2(npoints, &p1, npoints, &w1));
2003: for (PetscInt i = 0; i < totpoints * Nc; ++i) w[i] = 1.;
2004: for (PetscInt i = 0, totprev = 1, totrem = totpoints / npoints; i < dim; ++i) {
2005: PetscReal mul;
2007: mul = PetscPowReal(2., -i);
2008: PetscCall(PetscDTGaussJacobiQuadrature(npoints, -1., 1., i, 0.0, p1, w1));
2009: for (PetscInt pt = 0, l = 0; l < totprev; l++) {
2010: for (PetscInt j = 0; j < npoints; j++) {
2011: for (PetscInt m = 0; m < totrem; m++, pt++) {
2012: for (PetscInt k = 0; k < i; k++) x[pt * dim + k] = (x[pt * dim + k] + 1.) * (1. - p1[j]) * 0.5 - 1.;
2013: x[pt * dim + i] = p1[j];
2014: for (PetscInt c = 0; c < Nc; c++) w[pt * Nc + c] *= mul * w1[j];
2015: }
2016: }
2017: }
2018: totprev *= npoints;
2019: totrem /= npoints;
2020: }
2021: PetscCall(PetscFree2(p1, w1));
2022: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
2023: PetscCall(PetscQuadratureSetCellType(*q, ct));
2024: PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1));
2025: PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w));
2026: PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "StroudConical"));
2027: PetscFunctionReturn(PETSC_SUCCESS);
2028: }
2030: static PetscBool MinSymTriQuadCite = PETSC_FALSE;
2031: const char MinSymTriQuadCitation[] = "@article{WitherdenVincent2015,\n"
2032: " title = {On the identification of symmetric quadrature rules for finite element methods},\n"
2033: " journal = {Computers & Mathematics with Applications},\n"
2034: " volume = {69},\n"
2035: " number = {10},\n"
2036: " pages = {1232-1241},\n"
2037: " year = {2015},\n"
2038: " issn = {0898-1221},\n"
2039: " doi = {10.1016/j.camwa.2015.03.017},\n"
2040: " url = {https://www.sciencedirect.com/science/article/pii/S0898122115001224},\n"
2041: " author = {F.D. Witherden and P.E. Vincent},\n"
2042: "}\n";
2044: #include "petscdttriquadrules.h"
2046: static PetscBool MinSymTetQuadCite = PETSC_FALSE;
2047: const char MinSymTetQuadCitation[] = "@article{JaskowiecSukumar2021\n"
2048: " author = {Jaskowiec, Jan and Sukumar, N.},\n"
2049: " title = {High-order symmetric cubature rules for tetrahedra and pyramids},\n"
2050: " journal = {International Journal for Numerical Methods in Engineering},\n"
2051: " volume = {122},\n"
2052: " number = {1},\n"
2053: " pages = {148-171},\n"
2054: " doi = {10.1002/nme.6528},\n"
2055: " url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6528},\n"
2056: " eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6528},\n"
2057: " year = {2021}\n"
2058: "}\n";
2060: #include "petscdttetquadrules.h"
2062: static PetscBool DiagSymTriQuadCite = PETSC_FALSE;
2063: const char DiagSymTriQuadCitation[] = "@article{KongMulderVeldhuizen1999,\n"
2064: " title = {Higher-order triangular and tetrahedral finite elements with mass lumping for solving the wave equation},\n"
2065: " journal = {Journal of Engineering Mathematics},\n"
2066: " volume = {35},\n"
2067: " number = {4},\n"
2068: " pages = {405--426},\n"
2069: " year = {1999},\n"
2070: " doi = {10.1023/A:1004420829610},\n"
2071: " url = {https://link.springer.com/article/10.1023/A:1004420829610},\n"
2072: " author = {MJS Chin-Joe-Kong and Wim A Mulder and Marinus Van Veldhuizen},\n"
2073: "}\n";
2075: #include "petscdttridiagquadrules.h"
2077: // https://en.wikipedia.org/wiki/Partition_(number_theory)
2078: static PetscErrorCode PetscDTPartitionNumber(PetscInt n, PetscInt *p)
2079: {
2080: // sequence A000041 in the OEIS
2081: 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};
2082: PetscInt tabulated_max = PETSC_STATIC_ARRAY_LENGTH(partition) - 1;
2084: PetscFunctionBegin;
2085: PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Partition number not defined for negative number %" PetscInt_FMT, n);
2086: // not implementing the pentagonal number recurrence, we don't need partition numbers for n that high
2087: 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);
2088: *p = partition[n];
2089: PetscFunctionReturn(PETSC_SUCCESS);
2090: }
2092: /*@
2093: PetscDTSimplexQuadrature - Create a quadrature rule for a simplex that exactly integrates polynomials up to a given degree.
2095: Not Collective
2097: Input Parameters:
2098: + dim - The spatial dimension of the simplex (1 = segment, 2 = triangle, 3 = tetrahedron)
2099: . degree - The largest polynomial degree that is required to be integrated exactly
2100: - type - `PetscDTSimplexQuadratureType` indicating the type of quadrature rule
2102: Output Parameter:
2103: . quad - A `PetscQuadrature` object for integration over the biunit simplex
2104: (defined by the bounds $x_i >= -1$ and $\sum_i x_i <= 2 - d$) that is exact for
2105: polynomials up to the given degree
2107: Level: intermediate
2109: .seealso: `PetscDTSimplexQuadratureType`, `PetscDTGaussQuadrature()`, `PetscDTStroudCononicalQuadrature()`, `PetscQuadrature`
2110: @*/
2111: PetscErrorCode PetscDTSimplexQuadrature(PetscInt dim, PetscInt degree, PetscDTSimplexQuadratureType type, PetscQuadrature *quad)
2112: {
2113: PetscDTSimplexQuadratureType orig_type = type;
2115: PetscFunctionBegin;
2116: PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative dimension %" PetscInt_FMT, dim);
2117: PetscCheck(degree >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %" PetscInt_FMT, degree);
2118: if (type == PETSCDTSIMPLEXQUAD_DEFAULT) type = PETSCDTSIMPLEXQUAD_MINSYM;
2119: if (type == PETSCDTSIMPLEXQUAD_CONIC || dim < 2) {
2120: PetscInt points_per_dim = (degree + 2) / 2; // ceil((degree + 1) / 2);
2121: PetscCall(PetscDTStroudConicalQuadrature(dim, 1, points_per_dim, -1, 1, quad));
2122: } else {
2123: DMPolytopeType ct;
2124: PetscInt n = dim + 1;
2125: PetscInt fact = 1;
2126: PetscInt *part, *perm;
2127: PetscInt p = 0;
2128: PetscInt max_degree;
2129: const PetscInt *nodes_per_type = NULL;
2130: const PetscInt *all_num_full_nodes = NULL;
2131: const PetscReal **weights_list = NULL;
2132: const PetscReal **compact_nodes_list = NULL;
2133: const char *citation = NULL;
2134: PetscBool *cited = NULL;
2136: switch (dim) {
2137: case 0:
2138: ct = DM_POLYTOPE_POINT;
2139: break;
2140: case 1:
2141: ct = DM_POLYTOPE_SEGMENT;
2142: break;
2143: case 2:
2144: ct = DM_POLYTOPE_TRIANGLE;
2145: break;
2146: case 3:
2147: ct = DM_POLYTOPE_TETRAHEDRON;
2148: break;
2149: default:
2150: ct = DM_POLYTOPE_UNKNOWN;
2151: }
2152: if (type == PETSCDTSIMPLEXQUAD_MINSYM) {
2153: switch (dim) {
2154: case 2:
2155: cited = &MinSymTriQuadCite;
2156: citation = MinSymTriQuadCitation;
2157: max_degree = PetscDTWVTriQuad_max_degree;
2158: nodes_per_type = PetscDTWVTriQuad_num_orbits;
2159: all_num_full_nodes = PetscDTWVTriQuad_num_nodes;
2160: weights_list = PetscDTWVTriQuad_weights;
2161: compact_nodes_list = PetscDTWVTriQuad_orbits;
2162: break;
2163: case 3:
2164: cited = &MinSymTetQuadCite;
2165: citation = MinSymTetQuadCitation;
2166: max_degree = PetscDTJSTetQuad_max_degree;
2167: nodes_per_type = PetscDTJSTetQuad_num_orbits;
2168: all_num_full_nodes = PetscDTJSTetQuad_num_nodes;
2169: weights_list = PetscDTJSTetQuad_weights;
2170: compact_nodes_list = PetscDTJSTetQuad_orbits;
2171: break;
2172: default:
2173: max_degree = -1;
2174: break;
2175: }
2176: } else {
2177: switch (dim) {
2178: case 2:
2179: cited = &DiagSymTriQuadCite;
2180: citation = DiagSymTriQuadCitation;
2181: max_degree = PetscDTKMVTriQuad_max_degree;
2182: nodes_per_type = PetscDTKMVTriQuad_num_orbits;
2183: all_num_full_nodes = PetscDTKMVTriQuad_num_nodes;
2184: weights_list = PetscDTKMVTriQuad_weights;
2185: compact_nodes_list = PetscDTKMVTriQuad_orbits;
2186: break;
2187: default:
2188: max_degree = -1;
2189: break;
2190: }
2191: }
2193: if (degree > max_degree) {
2194: PetscCheck(orig_type == PETSCDTSIMPLEXQUAD_DEFAULT, PETSC_COMM_SELF, PETSC_ERR_SUP, "%s symmetric quadrature for dim %" PetscInt_FMT ", degree %" PetscInt_FMT " unsupported", orig_type == PETSCDTSIMPLEXQUAD_MINSYM ? "Minimal" : "Diagonal", dim, degree);
2195: // fall back to conic
2196: PetscCall(PetscDTSimplexQuadrature(dim, degree, PETSCDTSIMPLEXQUAD_CONIC, quad));
2197: PetscFunctionReturn(PETSC_SUCCESS);
2198: }
2200: PetscCall(PetscCitationsRegister(citation, cited));
2202: PetscCall(PetscDTPartitionNumber(n, &p));
2203: for (PetscInt d = 2; d <= n; d++) fact *= d;
2205: PetscInt num_full_nodes = all_num_full_nodes[degree];
2206: const PetscReal *all_compact_nodes = compact_nodes_list[degree];
2207: const PetscReal *all_compact_weights = weights_list[degree];
2208: nodes_per_type = &nodes_per_type[p * degree];
2210: PetscReal *points;
2211: PetscReal *counts;
2212: PetscReal *weights;
2213: PetscReal *bary_to_biunit; // row-major transformation of barycentric coordinate to biunit
2214: PetscQuadrature q;
2216: // compute the transformation
2217: PetscCall(PetscMalloc1(n * dim, &bary_to_biunit));
2218: for (PetscInt d = 0; d < dim; d++) {
2219: for (PetscInt b = 0; b < n; b++) bary_to_biunit[d * n + b] = (d == b) ? 1.0 : -1.0;
2220: }
2222: PetscCall(PetscMalloc3(n, &part, n, &perm, n, &counts));
2223: PetscCall(PetscCalloc1(num_full_nodes * dim, &points));
2224: PetscCall(PetscMalloc1(num_full_nodes, &weights));
2226: // (0, 0, ...) is the first partition lexicographically
2227: PetscCall(PetscArrayzero(part, n));
2228: PetscCall(PetscArrayzero(counts, n));
2229: counts[0] = n;
2231: // for each partition
2232: for (PetscInt s = 0, node_offset = 0; s < p; s++) {
2233: PetscInt num_compact_coords = part[n - 1] + 1;
2235: const PetscReal *compact_nodes = all_compact_nodes;
2236: const PetscReal *compact_weights = all_compact_weights;
2237: all_compact_nodes += num_compact_coords * nodes_per_type[s];
2238: all_compact_weights += nodes_per_type[s];
2240: // for every permutation of the vertices
2241: for (PetscInt f = 0; f < fact; f++) {
2242: PetscCall(PetscDTEnumPerm(n, f, perm, NULL));
2244: // check if it is a valid permutation
2245: PetscInt digit;
2246: for (digit = 1; digit < n; digit++) {
2247: // skip permutations that would duplicate a node because it has a smaller symmetry group
2248: if (part[digit - 1] == part[digit] && perm[digit - 1] > perm[digit]) break;
2249: }
2250: if (digit < n) continue;
2252: // create full nodes from this permutation of the compact nodes
2253: PetscReal *full_nodes = &points[node_offset * dim];
2254: PetscReal *full_weights = &weights[node_offset];
2256: PetscCall(PetscArraycpy(full_weights, compact_weights, nodes_per_type[s]));
2257: for (PetscInt b = 0; b < n; b++) {
2258: for (PetscInt d = 0; d < dim; d++) {
2259: 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]];
2260: }
2261: }
2262: node_offset += nodes_per_type[s];
2263: }
2265: if (s < p - 1) { // Generate the next partition
2266: /* A partition is described by the number of coordinates that are in
2267: * each set of duplicates (counts) and redundantly by mapping each
2268: * index to its set of duplicates (part)
2269: *
2270: * Counts should always be in nonincreasing order
2271: *
2272: * We want to generate the partitions lexically by part, which means
2273: * finding the last index where count > 1 and reducing by 1.
2274: *
2275: * For the new counts beyond that index, we eagerly assign the remaining
2276: * capacity of the partition to smaller indices (ensures lexical ordering),
2277: * while respecting the nonincreasing invariant of the counts
2278: */
2279: PetscInt last_digit = part[n - 1];
2280: PetscInt last_digit_with_extra = last_digit;
2281: while (counts[last_digit_with_extra] == 1) last_digit_with_extra--;
2282: PetscInt limit = --counts[last_digit_with_extra];
2283: PetscInt total_to_distribute = last_digit - last_digit_with_extra + 1;
2284: for (PetscInt digit = last_digit_with_extra + 1; digit < n; digit++) {
2285: counts[digit] = PetscMin(limit, total_to_distribute);
2286: total_to_distribute -= PetscMin(limit, total_to_distribute);
2287: }
2288: for (PetscInt digit = 0, offset = 0; digit < n; digit++) {
2289: PetscInt count = counts[digit];
2290: for (PetscInt c = 0; c < count; c++) part[offset++] = digit;
2291: }
2292: }
2293: PetscCheck(node_offset <= num_full_nodes, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node offset %" PetscInt_FMT " > %" PetscInt_FMT " number of nodes", node_offset, num_full_nodes);
2294: }
2295: PetscCall(PetscFree3(part, perm, counts));
2296: PetscCall(PetscFree(bary_to_biunit));
2297: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &q));
2298: PetscCall(PetscQuadratureSetCellType(q, ct));
2299: PetscCall(PetscQuadratureSetOrder(q, degree));
2300: PetscCall(PetscQuadratureSetData(q, dim, 1, num_full_nodes, points, weights));
2301: *quad = q;
2302: }
2303: PetscFunctionReturn(PETSC_SUCCESS);
2304: }
2306: /*@
2307: PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell
2309: Not Collective
2311: Input Parameters:
2312: + dim - The cell dimension
2313: . level - The number of points in one dimension, $2^l$
2314: . a - left end of interval (often-1)
2315: - b - right end of interval (often +1)
2317: Output Parameter:
2318: . q - A `PetscQuadrature` object
2320: Level: intermediate
2322: .seealso: `PetscDTGaussTensorQuadrature()`, `PetscQuadrature`
2323: @*/
2324: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
2325: {
2326: DMPolytopeType ct;
2327: const PetscInt p = 16; /* Digits of precision in the evaluation */
2328: const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */
2329: const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */
2330: const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */
2331: PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */
2332: PetscReal wk = 0.5 * PETSC_PI; /* Quadrature weight at x_k */
2333: PetscReal *x, *w;
2334: PetscInt K, k, npoints;
2336: PetscFunctionBegin;
2337: PetscCheck(dim <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not yet implemented", dim);
2338: PetscCheck(level, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
2339: switch (dim) {
2340: case 0:
2341: ct = DM_POLYTOPE_POINT;
2342: break;
2343: case 1:
2344: ct = DM_POLYTOPE_SEGMENT;
2345: break;
2346: case 2:
2347: ct = DM_POLYTOPE_QUADRILATERAL;
2348: break;
2349: case 3:
2350: ct = DM_POLYTOPE_HEXAHEDRON;
2351: break;
2352: default:
2353: ct = DM_POLYTOPE_UNKNOWN;
2354: }
2355: /* Find K such that the weights are < 32 digits of precision */
2356: 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)));
2357: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
2358: PetscCall(PetscQuadratureSetCellType(*q, ct));
2359: PetscCall(PetscQuadratureSetOrder(*q, 2 * K + 1));
2360: npoints = 2 * K - 1;
2361: PetscCall(PetscMalloc1(npoints * dim, &x));
2362: PetscCall(PetscMalloc1(npoints, &w));
2363: /* Center term */
2364: x[0] = beta;
2365: w[0] = 0.5 * alpha * h * PETSC_PI;
2366: for (k = 1; k < K; ++k) {
2367: wk = 0.5 * alpha * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2368: xk = PetscTanhReal(0.5 * PETSC_PI * PetscSinhReal(k * h));
2369: x[2 * k - 1] = -alpha * xk + beta;
2370: w[2 * k - 1] = wk;
2371: x[2 * k + 0] = alpha * xk + beta;
2372: w[2 * k + 0] = wk;
2373: }
2374: PetscCall(PetscQuadratureSetData(*q, dim, 1, npoints, x, w));
2375: PetscFunctionReturn(PETSC_SUCCESS);
2376: }
2378: /*@
2379: PetscDTTanhSinhIntegrate - Approximate $\int_a^b f(x)\,dx$ to a requested precision using adaptive tanh-sinh (double-exponential) quadrature
2381: Not Collective; No Fortran Support
2383: Input Parameters:
2384: + func - the integrand callback (`func(x, ctx, &value)` evaluates the integrand at point `x`)
2385: . a - lower limit of integration
2386: . b - upper limit of integration
2387: . digits - target number of correct decimal digits
2388: - ctx - optional application context passed to `func`
2390: Output Parameter:
2391: . sol - the approximate value of the integral
2393: Level: developer
2395: Notes:
2396: Doubles the number of quadrature points at each refinement level until the change in the
2397: integral falls below the requested tolerance. Suitable for smooth integrands and integrands
2398: with endpoint singularities.
2400: For arbitrary-precision arithmetic via MPFR, see `PetscDTTanhSinhIntegrateMPFR()`.
2402: .seealso: `PetscDTTanhSinhIntegrateMPFR()`, `PetscDTGaussQuadrature()`
2403: @*/
2404: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(const PetscReal[], PetscCtx, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscCtx ctx, PetscReal *sol)
2405: {
2406: const PetscInt p = 16; /* Digits of precision in the evaluation */
2407: const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */
2408: const PetscReal beta = (b + a) / 2.; /* Center of the integration interval */
2409: PetscReal h = 1.0; /* Step size, length between x_k */
2410: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
2411: PetscReal osum = 0.0; /* Integral on last level */
2412: PetscReal psum = 0.0; /* Integral on the level before the last level */
2413: PetscReal sum; /* Integral on current level */
2414: PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
2415: PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
2416: PetscReal wk; /* Quadrature weight at x_k */
2417: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
2418: PetscInt d; /* Digits of precision in the integral */
2420: PetscFunctionBegin;
2421: PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
2422: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2423: /* Center term */
2424: func(&beta, ctx, &lval);
2425: sum = 0.5 * alpha * PETSC_PI * lval;
2426: /* */
2427: do {
2428: PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
2429: PetscInt k = 1;
2431: ++l;
2432: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */
2433: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
2434: psum = osum;
2435: osum = sum;
2436: h *= 0.5;
2437: sum *= 0.5;
2438: do {
2439: wk = 0.5 * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2440: yk = 1.0 / (PetscExpReal(0.5 * PETSC_PI * PetscSinhReal(k * h)) * PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2441: lx = -alpha * (1.0 - yk) + beta;
2442: rx = alpha * (1.0 - yk) + beta;
2443: func(&lx, ctx, &lval);
2444: func(&rx, ctx, &rval);
2445: lterm = alpha * wk * lval;
2446: maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
2447: sum += lterm;
2448: rterm = alpha * wk * rval;
2449: maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
2450: sum += rterm;
2451: ++k;
2452: /* Only need to evaluate every other point on refined levels */
2453: if (l != 1) ++k;
2454: } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */
2456: d1 = PetscLog10Real(PetscAbsReal(sum - osum));
2457: d2 = PetscLog10Real(PetscAbsReal(sum - psum));
2458: d3 = PetscLog10Real(maxTerm) - p;
2459: if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
2460: else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
2461: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4)));
2462: } while (d < digits && l < 12);
2463: *sol = sum;
2464: PetscCall(PetscFPTrapPop());
2465: PetscFunctionReturn(PETSC_SUCCESS);
2466: }
2468: /*@
2469: PetscDTTanhSinhIntegrateMPFR - High-precision version of `PetscDTTanhSinhIntegrate()` that uses the MPFR arbitrary-precision library to evaluate the quadrature
2471: Not Collective; No Fortran Support
2473: Input Parameters:
2474: + func - the integrand callback (`func(x, ctx, &value)` evaluates the integrand at point `x`)
2475: . a - lower limit of integration
2476: . b - upper limit of integration
2477: . digits - target number of correct decimal digits (also drives the working MPFR precision)
2478: - ctx - optional application context passed to `func`
2480: Output Parameter:
2481: . sol - the approximate value of the integral
2483: Level: developer
2485: Note:
2486: Requires PETSc to be configured with `--with-mpfr`; otherwise an error is raised.
2488: .seealso: `PetscDTTanhSinhIntegrate()`, `PetscDTGaussQuadrature()`
2489: @*/
2490: #if PetscDefined(HAVE_MPFR)
2491: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], PetscCtx, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscCtx ctx, PetscReal *sol)
2492: {
2493: const PetscInt safetyFactor = 2; /* Calculate abscissa until 2*p digits */
2494: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
2495: mpfr_t alpha; /* Half-width of the integration interval */
2496: mpfr_t beta; /* Center of the integration interval */
2497: mpfr_t h; /* Step size, length between x_k */
2498: mpfr_t osum; /* Integral on last level */
2499: mpfr_t psum; /* Integral on the level before the last level */
2500: mpfr_t sum; /* Integral on current level */
2501: mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
2502: mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
2503: mpfr_t wk; /* Quadrature weight at x_k */
2504: PetscReal lval, rval, rtmp; /* Terms in the quadature sum to the left and right of 0 */
2505: PetscInt d; /* Digits of precision in the integral */
2506: mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;
2508: PetscFunctionBegin;
2509: PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
2510: /* Create high precision storage */
2511: 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);
2512: /* Initialization */
2513: mpfr_set_d(alpha, 0.5 * (b - a), MPFR_RNDN);
2514: mpfr_set_d(beta, 0.5 * (b + a), MPFR_RNDN);
2515: mpfr_set_d(osum, 0.0, MPFR_RNDN);
2516: mpfr_set_d(psum, 0.0, MPFR_RNDN);
2517: mpfr_set_d(h, 1.0, MPFR_RNDN);
2518: mpfr_const_pi(pi2, MPFR_RNDN);
2519: mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
2520: /* Center term */
2521: rtmp = 0.5 * (b + a);
2522: func(&rtmp, ctx, &lval);
2523: mpfr_set(sum, pi2, MPFR_RNDN);
2524: mpfr_mul(sum, sum, alpha, MPFR_RNDN);
2525: mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
2526: /* */
2527: do {
2528: PetscReal d1, d2, d3, d4;
2529: PetscInt k = 1;
2531: ++l;
2532: mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
2533: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */
2534: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
2535: mpfr_set(psum, osum, MPFR_RNDN);
2536: mpfr_set(osum, sum, MPFR_RNDN);
2537: mpfr_mul_d(h, h, 0.5, MPFR_RNDN);
2538: mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
2539: do {
2540: mpfr_set_si(kh, k, MPFR_RNDN);
2541: mpfr_mul(kh, kh, h, MPFR_RNDN);
2542: /* Weight */
2543: mpfr_set(wk, h, MPFR_RNDN);
2544: mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
2545: mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
2546: mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
2547: mpfr_cosh(tmp, msinh, MPFR_RNDN);
2548: mpfr_sqr(tmp, tmp, MPFR_RNDN);
2549: mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
2550: mpfr_div(wk, wk, tmp, MPFR_RNDN);
2551: /* Abscissa */
2552: mpfr_set_d(yk, 1.0, MPFR_RNDZ);
2553: mpfr_cosh(tmp, msinh, MPFR_RNDN);
2554: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
2555: mpfr_exp(tmp, msinh, MPFR_RNDN);
2556: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
2557: /* Quadrature points */
2558: mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
2559: mpfr_mul(lx, lx, alpha, MPFR_RNDU);
2560: mpfr_add(lx, lx, beta, MPFR_RNDU);
2561: mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
2562: mpfr_mul(rx, rx, alpha, MPFR_RNDD);
2563: mpfr_add(rx, rx, beta, MPFR_RNDD);
2564: /* Evaluation */
2565: rtmp = mpfr_get_d(lx, MPFR_RNDU);
2566: func(&rtmp, ctx, &lval);
2567: rtmp = mpfr_get_d(rx, MPFR_RNDD);
2568: func(&rtmp, ctx, &rval);
2569: /* Update */
2570: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
2571: mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
2572: mpfr_add(sum, sum, tmp, MPFR_RNDN);
2573: mpfr_abs(tmp, tmp, MPFR_RNDN);
2574: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
2575: mpfr_set(curTerm, tmp, MPFR_RNDN);
2576: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
2577: mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
2578: mpfr_add(sum, sum, tmp, MPFR_RNDN);
2579: mpfr_abs(tmp, tmp, MPFR_RNDN);
2580: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
2581: mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
2582: ++k;
2583: /* Only need to evaluate every other point on refined levels */
2584: if (l != 1) ++k;
2585: mpfr_log10(tmp, wk, MPFR_RNDN);
2586: mpfr_abs(tmp, tmp, MPFR_RNDN);
2587: } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor * digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
2588: mpfr_sub(tmp, sum, osum, MPFR_RNDN);
2589: mpfr_abs(tmp, tmp, MPFR_RNDN);
2590: mpfr_log10(tmp, tmp, MPFR_RNDN);
2591: d1 = mpfr_get_d(tmp, MPFR_RNDN);
2592: mpfr_sub(tmp, sum, psum, MPFR_RNDN);
2593: mpfr_abs(tmp, tmp, MPFR_RNDN);
2594: mpfr_log10(tmp, tmp, MPFR_RNDN);
2595: d2 = mpfr_get_d(tmp, MPFR_RNDN);
2596: mpfr_log10(tmp, maxTerm, MPFR_RNDN);
2597: d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
2598: mpfr_log10(tmp, curTerm, MPFR_RNDN);
2599: d4 = mpfr_get_d(tmp, MPFR_RNDN);
2600: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4)));
2601: } while (d < digits && l < 8);
2602: *sol = mpfr_get_d(sum, MPFR_RNDN);
2603: /* Cleanup */
2604: mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
2605: PetscFunctionReturn(PETSC_SUCCESS);
2606: }
2607: #else
2609: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscCtx ctx, PetscReal *sol)
2610: {
2611: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
2612: }
2613: #endif
2615: /*@
2616: PetscDTTensorQuadratureCreate - create the tensor product quadrature from two lower-dimensional quadratures
2618: Not Collective
2620: Input Parameters:
2621: + q1 - The first quadrature
2622: - q2 - The second quadrature
2624: Output Parameter:
2625: . q - A `PetscQuadrature` object
2627: Level: intermediate
2629: .seealso: `PetscQuadrature`, `PetscDTGaussTensorQuadrature()`
2630: @*/
2631: PetscErrorCode PetscDTTensorQuadratureCreate(PetscQuadrature q1, PetscQuadrature q2, PetscQuadrature *q)
2632: {
2633: DMPolytopeType ct1, ct2, ct;
2634: const PetscReal *x1, *w1, *x2, *w2;
2635: PetscReal *x, *w;
2636: PetscInt dim1, Nc1, Np1, order1, qa, d1;
2637: PetscInt dim2, Nc2, Np2, order2, qb, d2;
2638: PetscInt dim, Nc, Np, order, qc, d;
2640: PetscFunctionBegin;
2643: PetscAssertPointer(q, 3);
2645: PetscCall(PetscQuadratureGetOrder(q1, &order1));
2646: PetscCall(PetscQuadratureGetOrder(q2, &order2));
2647: PetscCheck(order1 == order2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Order1 %" PetscInt_FMT " != %" PetscInt_FMT " Order2", order1, order2);
2648: PetscCall(PetscQuadratureGetData(q1, &dim1, &Nc1, &Np1, &x1, &w1));
2649: PetscCall(PetscQuadratureGetCellType(q1, &ct1));
2650: PetscCall(PetscQuadratureGetData(q2, &dim2, &Nc2, &Np2, &x2, &w2));
2651: PetscCall(PetscQuadratureGetCellType(q2, &ct2));
2652: PetscCheck(Nc1 == Nc2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "NumComp1 %" PetscInt_FMT " != %" PetscInt_FMT " NumComp2", Nc1, Nc2);
2654: switch (ct1) {
2655: case DM_POLYTOPE_POINT:
2656: ct = ct2;
2657: break;
2658: case DM_POLYTOPE_SEGMENT:
2659: switch (ct2) {
2660: case DM_POLYTOPE_POINT:
2661: ct = ct1;
2662: break;
2663: case DM_POLYTOPE_SEGMENT:
2664: ct = DM_POLYTOPE_QUADRILATERAL;
2665: break;
2666: case DM_POLYTOPE_TRIANGLE:
2667: ct = DM_POLYTOPE_TRI_PRISM;
2668: break;
2669: case DM_POLYTOPE_QUADRILATERAL:
2670: ct = DM_POLYTOPE_HEXAHEDRON;
2671: break;
2672: case DM_POLYTOPE_TETRAHEDRON:
2673: ct = DM_POLYTOPE_UNKNOWN;
2674: break;
2675: case DM_POLYTOPE_HEXAHEDRON:
2676: ct = DM_POLYTOPE_UNKNOWN;
2677: break;
2678: default:
2679: ct = DM_POLYTOPE_UNKNOWN;
2680: }
2681: break;
2682: case DM_POLYTOPE_TRIANGLE:
2683: switch (ct2) {
2684: case DM_POLYTOPE_POINT:
2685: ct = ct1;
2686: break;
2687: case DM_POLYTOPE_SEGMENT:
2688: ct = DM_POLYTOPE_TRI_PRISM;
2689: break;
2690: case DM_POLYTOPE_TRIANGLE:
2691: ct = DM_POLYTOPE_UNKNOWN;
2692: break;
2693: case DM_POLYTOPE_QUADRILATERAL:
2694: ct = DM_POLYTOPE_UNKNOWN;
2695: break;
2696: case DM_POLYTOPE_TETRAHEDRON:
2697: ct = DM_POLYTOPE_UNKNOWN;
2698: break;
2699: case DM_POLYTOPE_HEXAHEDRON:
2700: ct = DM_POLYTOPE_UNKNOWN;
2701: break;
2702: default:
2703: ct = DM_POLYTOPE_UNKNOWN;
2704: }
2705: break;
2706: case DM_POLYTOPE_QUADRILATERAL:
2707: switch (ct2) {
2708: case DM_POLYTOPE_POINT:
2709: ct = ct1;
2710: break;
2711: case DM_POLYTOPE_SEGMENT:
2712: ct = DM_POLYTOPE_HEXAHEDRON;
2713: break;
2714: case DM_POLYTOPE_TRIANGLE:
2715: ct = DM_POLYTOPE_UNKNOWN;
2716: break;
2717: case DM_POLYTOPE_QUADRILATERAL:
2718: ct = DM_POLYTOPE_UNKNOWN;
2719: break;
2720: case DM_POLYTOPE_TETRAHEDRON:
2721: ct = DM_POLYTOPE_UNKNOWN;
2722: break;
2723: case DM_POLYTOPE_HEXAHEDRON:
2724: ct = DM_POLYTOPE_UNKNOWN;
2725: break;
2726: default:
2727: ct = DM_POLYTOPE_UNKNOWN;
2728: }
2729: break;
2730: case DM_POLYTOPE_TETRAHEDRON:
2731: switch (ct2) {
2732: case DM_POLYTOPE_POINT:
2733: ct = ct1;
2734: break;
2735: case DM_POLYTOPE_SEGMENT:
2736: ct = DM_POLYTOPE_UNKNOWN;
2737: break;
2738: case DM_POLYTOPE_TRIANGLE:
2739: ct = DM_POLYTOPE_UNKNOWN;
2740: break;
2741: case DM_POLYTOPE_QUADRILATERAL:
2742: ct = DM_POLYTOPE_UNKNOWN;
2743: break;
2744: case DM_POLYTOPE_TETRAHEDRON:
2745: ct = DM_POLYTOPE_UNKNOWN;
2746: break;
2747: case DM_POLYTOPE_HEXAHEDRON:
2748: ct = DM_POLYTOPE_UNKNOWN;
2749: break;
2750: default:
2751: ct = DM_POLYTOPE_UNKNOWN;
2752: }
2753: break;
2754: case DM_POLYTOPE_HEXAHEDRON:
2755: switch (ct2) {
2756: case DM_POLYTOPE_POINT:
2757: ct = ct1;
2758: break;
2759: case DM_POLYTOPE_SEGMENT:
2760: ct = DM_POLYTOPE_UNKNOWN;
2761: break;
2762: case DM_POLYTOPE_TRIANGLE:
2763: ct = DM_POLYTOPE_UNKNOWN;
2764: break;
2765: case DM_POLYTOPE_QUADRILATERAL:
2766: ct = DM_POLYTOPE_UNKNOWN;
2767: break;
2768: case DM_POLYTOPE_TETRAHEDRON:
2769: ct = DM_POLYTOPE_UNKNOWN;
2770: break;
2771: case DM_POLYTOPE_HEXAHEDRON:
2772: ct = DM_POLYTOPE_UNKNOWN;
2773: break;
2774: default:
2775: ct = DM_POLYTOPE_UNKNOWN;
2776: }
2777: break;
2778: default:
2779: ct = DM_POLYTOPE_UNKNOWN;
2780: }
2781: dim = dim1 + dim2;
2782: Nc = Nc1;
2783: Np = Np1 * Np2;
2784: order = order1;
2785: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
2786: PetscCall(PetscQuadratureSetCellType(*q, ct));
2787: PetscCall(PetscQuadratureSetOrder(*q, order));
2788: PetscCall(PetscMalloc1(Np * dim, &x));
2789: PetscCall(PetscMalloc1(Np, &w));
2790: for (qa = 0, qc = 0; qa < Np1; ++qa) {
2791: for (qb = 0; qb < Np2; ++qb, ++qc) {
2792: for (d1 = 0, d = 0; d1 < dim1; ++d1, ++d) x[qc * dim + d] = x1[qa * dim1 + d1];
2793: for (d2 = 0; d2 < dim2; ++d2, ++d) x[qc * dim + d] = x2[qb * dim2 + d2];
2794: w[qc] = w1[qa] * w2[qb];
2795: }
2796: }
2797: PetscCall(PetscQuadratureSetData(*q, dim, Nc, Np, x, w));
2798: PetscFunctionReturn(PETSC_SUCCESS);
2799: }
2801: /* Overwrites A. Can only handle full-rank problems with m>=n
2802: A in column-major format
2803: Ainv in row-major format
2804: tau has length m
2805: worksize must be >= max(1,n)
2806: */
2807: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m, PetscInt mstride, PetscInt n, PetscReal *A_in, PetscReal *Ainv_out, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2808: {
2809: PetscBLASInt M, N, K, lda, ldb, ldwork;
2810: PetscScalar *A, *Ainv, *R, *Q, Alpha;
2812: PetscFunctionBegin;
2813: #if PetscDefined(USE_COMPLEX)
2814: {
2815: PetscInt i, j;
2816: PetscCall(PetscMalloc2(m * n, &A, m * n, &Ainv));
2817: for (j = 0; j < n; j++) {
2818: for (i = 0; i < m; i++) A[i + m * j] = A_in[i + mstride * j];
2819: }
2820: mstride = m;
2821: }
2822: #else
2823: A = A_in;
2824: Ainv = Ainv_out;
2825: #endif
2827: PetscCall(PetscBLASIntCast(m, &M));
2828: PetscCall(PetscBLASIntCast(n, &N));
2829: PetscCall(PetscBLASIntCast(mstride, &lda));
2830: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2831: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2832: PetscCallLAPACKInfo("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
2833: PetscCall(PetscFPTrapPop());
2834: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2836: /* Extract an explicit representation of Q */
2837: Q = Ainv;
2838: PetscCall(PetscArraycpy(Q, A, mstride * n));
2839: K = N; /* full rank */
2840: PetscCallLAPACKInfo("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
2842: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2843: Alpha = 1.0;
2844: ldb = lda;
2845: PetscCallBLAS("BLAStrsm", BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb));
2846: /* Ainv is Q, overwritten with inverse */
2848: #if PetscDefined(USE_COMPLEX)
2849: {
2850: PetscInt i;
2851: for (i = 0; i < m * n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
2852: PetscCall(PetscFree2(A, Ainv));
2853: }
2854: #endif
2855: PetscFunctionReturn(PETSC_SUCCESS);
2856: }
2858: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
2859: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval, const PetscReal *x, PetscInt ndegree, const PetscInt *degrees, PetscBool Transpose, PetscReal *B)
2860: {
2861: PetscReal *Bv;
2862: PetscInt i, j;
2864: PetscFunctionBegin;
2865: PetscCall(PetscMalloc1((ninterval + 1) * ndegree, &Bv));
2866: /* Point evaluation of L_p on all the source vertices */
2867: PetscCall(PetscDTLegendreEval(ninterval + 1, x, ndegree, degrees, Bv, NULL, NULL));
2868: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
2869: for (i = 0; i < ninterval; i++) {
2870: for (j = 0; j < ndegree; j++) {
2871: if (Transpose) B[i + ninterval * j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j];
2872: else B[i * ndegree + j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j];
2873: }
2874: }
2875: PetscCall(PetscFree(Bv));
2876: PetscFunctionReturn(PETSC_SUCCESS);
2877: }
2879: /*@
2880: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
2882: Not Collective
2884: Input Parameters:
2885: + degree - degree of reconstruction polynomial
2886: . nsource - number of source intervals
2887: . sourcex - sorted coordinates of source cell boundaries (length `nsource`+1)
2888: . ntarget - number of target intervals
2889: - targetx - sorted coordinates of target cell boundaries (length `ntarget`+1)
2891: Output Parameter:
2892: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
2894: Level: advanced
2896: .seealso: `PetscDTLegendreEval()`
2897: @*/
2898: PetscErrorCode PetscDTReconstructPoly(PetscInt degree, PetscInt nsource, const PetscReal sourcex[], PetscInt ntarget, const PetscReal targetx[], PetscReal R[])
2899: {
2900: PetscInt i, j, k, *bdegrees, worksize;
2901: PetscReal xmin, xmax, center, hscale, *sourcey, *targety, *Bsource, *Bsinv, *Btarget;
2902: PetscScalar *tau, *work;
2904: PetscFunctionBegin;
2905: PetscAssertPointer(sourcex, 3);
2906: PetscAssertPointer(targetx, 5);
2907: PetscAssertPointer(R, 6);
2908: 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);
2909: if (PetscDefined(USE_DEBUG)) {
2910: 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]);
2911: 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]);
2912: }
2913: xmin = PetscMin(sourcex[0], targetx[0]);
2914: xmax = PetscMax(sourcex[nsource], targetx[ntarget]);
2915: center = (xmin + xmax) / 2;
2916: hscale = (xmax - xmin) / 2;
2917: worksize = nsource;
2918: PetscCall(PetscMalloc4(degree + 1, &bdegrees, nsource + 1, &sourcey, nsource * (degree + 1), &Bsource, worksize, &work));
2919: PetscCall(PetscMalloc4(nsource, &tau, nsource * (degree + 1), &Bsinv, ntarget + 1, &targety, ntarget * (degree + 1), &Btarget));
2920: for (i = 0; i <= nsource; i++) sourcey[i] = (sourcex[i] - center) / hscale;
2921: for (i = 0; i <= degree; i++) bdegrees[i] = i + 1;
2922: PetscCall(PetscDTLegendreIntegrate(nsource, sourcey, degree + 1, bdegrees, PETSC_TRUE, Bsource));
2923: PetscCall(PetscDTPseudoInverseQR(nsource, nsource, degree + 1, Bsource, Bsinv, tau, nsource, work));
2924: for (i = 0; i <= ntarget; i++) targety[i] = (targetx[i] - center) / hscale;
2925: PetscCall(PetscDTLegendreIntegrate(ntarget, targety, degree + 1, bdegrees, PETSC_FALSE, Btarget));
2926: for (i = 0; i < ntarget; i++) {
2927: PetscReal rowsum = 0;
2928: for (j = 0; j < nsource; j++) {
2929: PetscReal sum = 0;
2930: for (k = 0; k < degree + 1; k++) sum += Btarget[i * (degree + 1) + k] * Bsinv[k * nsource + j];
2931: R[i * nsource + j] = sum;
2932: rowsum += sum;
2933: }
2934: for (j = 0; j < nsource; j++) R[i * nsource + j] /= rowsum; /* normalize each row */
2935: }
2936: PetscCall(PetscFree4(bdegrees, sourcey, Bsource, work));
2937: PetscCall(PetscFree4(tau, Bsinv, targety, Btarget));
2938: PetscFunctionReturn(PETSC_SUCCESS);
2939: }
2941: /*@
2942: PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points
2944: Not Collective
2946: Input Parameters:
2947: + n - the number of GLL nodes
2948: . nodes - the GLL nodes
2949: . weights - the GLL weights
2950: - f - the function values at the nodes
2952: Output Parameter:
2953: . in - the value of the integral
2955: Level: beginner
2957: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`
2958: @*/
2959: PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n, PetscReal nodes[], PetscReal weights[], const PetscReal f[], PetscReal *in)
2960: {
2961: PetscInt i;
2963: PetscFunctionBegin;
2964: *in = 0.;
2965: for (i = 0; i < n; i++) *in += f[i] * f[i] * weights[i];
2966: PetscFunctionReturn(PETSC_SUCCESS);
2967: }
2969: /*@
2970: PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element
2972: Not Collective
2974: Input Parameters:
2975: + n - the number of GLL nodes
2976: . nodes - the GLL nodes, of length `n`
2977: - weights - the GLL weights, of length `n`
2979: Output Parameter:
2980: . AA - the stiffness element, of size `n` by `n`
2982: Level: beginner
2984: Notes:
2985: Destroy this with `PetscGaussLobattoLegendreElementLaplacianDestroy()`
2987: 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)
2989: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`
2990: @*/
2991: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA)
2992: {
2993: PetscReal **A;
2994: const PetscReal *gllnodes = nodes;
2995: const PetscInt p = n - 1;
2996: PetscReal z0, z1, z2 = -1, x, Lpj, Lpr;
2997: PetscInt i, j, nn, r;
2999: PetscFunctionBegin;
3000: PetscCall(PetscMalloc1(n, &A));
3001: PetscCall(PetscMalloc1(n * n, &A[0]));
3002: for (i = 1; i < n; i++) A[i] = A[i - 1] + n;
3004: for (j = 1; j < p; j++) {
3005: x = gllnodes[j];
3006: z0 = 1.;
3007: z1 = x;
3008: for (nn = 1; nn < p; nn++) {
3009: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
3010: z0 = z1;
3011: z1 = z2;
3012: }
3013: Lpj = z2;
3014: for (r = 1; r < p; r++) {
3015: if (r == j) {
3016: A[j][j] = 2. / (3. * (1. - gllnodes[j] * gllnodes[j]) * Lpj * Lpj);
3017: } else {
3018: x = gllnodes[r];
3019: z0 = 1.;
3020: z1 = x;
3021: for (nn = 1; nn < p; nn++) {
3022: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
3023: z0 = z1;
3024: z1 = z2;
3025: }
3026: Lpr = z2;
3027: A[r][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * Lpr * (gllnodes[j] - gllnodes[r]) * (gllnodes[j] - gllnodes[r]));
3028: }
3029: }
3030: }
3031: for (j = 1; j < p + 1; j++) {
3032: x = gllnodes[j];
3033: z0 = 1.;
3034: z1 = x;
3035: for (nn = 1; nn < p; nn++) {
3036: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
3037: z0 = z1;
3038: z1 = z2;
3039: }
3040: Lpj = z2;
3041: A[j][0] = 4. * PetscPowRealInt(-1., p) / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. + gllnodes[j]) * (1. + gllnodes[j]));
3042: A[0][j] = A[j][0];
3043: }
3044: for (j = 0; j < p; j++) {
3045: x = gllnodes[j];
3046: z0 = 1.;
3047: z1 = x;
3048: for (nn = 1; nn < p; nn++) {
3049: z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
3050: z0 = z1;
3051: z1 = z2;
3052: }
3053: Lpj = z2;
3055: A[p][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. - gllnodes[j]) * (1. - gllnodes[j]));
3056: A[j][p] = A[p][j];
3057: }
3058: A[0][0] = 0.5 + (((PetscReal)p) * (((PetscReal)p) + 1.) - 2.) / 6.;
3059: A[p][p] = A[0][0];
3060: *AA = A;
3061: PetscFunctionReturn(PETSC_SUCCESS);
3062: }
3064: /*@
3065: PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element created with `PetscGaussLobattoLegendreElementLaplacianCreate()`
3067: Not Collective
3069: Input Parameters:
3070: + n - the number of GLL nodes
3071: . nodes - the GLL nodes, ignored
3072: . weights - the GLL weightss, ignored
3073: - AA - the stiffness element from `PetscGaussLobattoLegendreElementLaplacianCreate()`
3075: Level: beginner
3077: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`
3078: @*/
3079: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA)
3080: {
3081: PetscFunctionBegin;
3082: PetscCall(PetscFree((*AA)[0]));
3083: PetscCall(PetscFree(*AA));
3084: *AA = NULL;
3085: PetscFunctionReturn(PETSC_SUCCESS);
3086: }
3088: /*@
3089: PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element
3091: Not Collective
3093: Input Parameters:
3094: + n - the number of GLL nodes
3095: . nodes - the GLL nodes, of length `n`
3096: - weights - the GLL weights, of length `n`
3098: Output Parameters:
3099: + AA - the stiffness element, of dimension `n` by `n`
3100: - AAT - the transpose of AA (pass in `NULL` if you do not need this array), of dimension `n` by `n`
3102: Level: beginner
3104: Notes:
3105: Destroy this with `PetscGaussLobattoLegendreElementGradientDestroy()`
3107: You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row-oriented
3109: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`, `PetscGaussLobattoLegendreElementGradientDestroy()`
3110: @*/
3111: PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA, PetscReal ***AAT)
3112: {
3113: PetscReal **A, **AT = NULL;
3114: const PetscReal *gllnodes = nodes;
3115: const PetscInt p = n - 1;
3116: PetscReal Li, Lj, d0;
3117: PetscInt i, j;
3119: PetscFunctionBegin;
3120: PetscCall(PetscMalloc1(n, &A));
3121: PetscCall(PetscMalloc1(n * n, &A[0]));
3122: for (i = 1; i < n; i++) A[i] = A[i - 1] + n;
3124: if (AAT) {
3125: PetscCall(PetscMalloc1(n, &AT));
3126: PetscCall(PetscMalloc1(n * n, &AT[0]));
3127: for (i = 1; i < n; i++) AT[i] = AT[i - 1] + n;
3128: }
3130: if (n == 1) A[0][0] = 0.;
3131: d0 = (PetscReal)p * ((PetscReal)p + 1.) / 4.;
3132: for (i = 0; i < n; i++) {
3133: for (j = 0; j < n; j++) {
3134: A[i][j] = 0.;
3135: PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li));
3136: PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj));
3137: if (i != j) A[i][j] = Li / (Lj * (gllnodes[i] - gllnodes[j]));
3138: if (j == i && i == 0) A[i][j] = -d0;
3139: if (j == i && i == p) A[i][j] = d0;
3140: if (AT) AT[j][i] = A[i][j];
3141: }
3142: }
3143: if (AAT) *AAT = AT;
3144: *AA = A;
3145: PetscFunctionReturn(PETSC_SUCCESS);
3146: }
3148: /*@
3149: PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with `PetscGaussLobattoLegendreElementGradientCreate()`
3151: Not Collective
3153: Input Parameters:
3154: + n - the number of GLL nodes
3155: . nodes - the GLL nodes, ignored
3156: . weights - the GLL weights, ignored
3157: . AA - the stiffness element obtained with `PetscGaussLobattoLegendreElementGradientCreate()`
3158: - AAT - the transpose of the element obtained with `PetscGaussLobattoLegendreElementGradientCreate()`
3160: Level: beginner
3162: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionCreate()`
3163: @*/
3164: PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA, PetscReal ***AAT)
3165: {
3166: PetscFunctionBegin;
3167: PetscCall(PetscFree((*AA)[0]));
3168: PetscCall(PetscFree(*AA));
3169: *AA = NULL;
3170: if (AAT) {
3171: PetscCall(PetscFree((*AAT)[0]));
3172: PetscCall(PetscFree(*AAT));
3173: *AAT = NULL;
3174: }
3175: PetscFunctionReturn(PETSC_SUCCESS);
3176: }
3178: /*@
3179: PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element
3181: Not Collective
3183: Input Parameters:
3184: + n - the number of GLL nodes
3185: . nodes - the GLL nodes, of length `n`
3186: - weights - the GLL weights, of length `n`
3188: Output Parameter:
3189: . AA - the stiffness element, of dimension `n` by `n`
3191: Level: beginner
3193: Notes:
3194: Destroy this with `PetscGaussLobattoLegendreElementAdvectionDestroy()`
3196: This is the same as the Gradient operator multiplied by the diagonal mass matrix
3198: You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row-oriented
3200: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionDestroy()`
3201: @*/
3202: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA)
3203: {
3204: PetscReal **D;
3205: const PetscReal *gllweights = weights;
3206: const PetscInt glln = n;
3207: PetscInt j;
3209: PetscFunctionBegin;
3210: PetscCall(PetscGaussLobattoLegendreElementGradientCreate(n, nodes, weights, &D, NULL));
3211: for (PetscInt i = 0; i < glln; i++) {
3212: for (j = 0; j < glln; j++) D[i][j] = gllweights[i] * D[i][j];
3213: }
3214: *AA = D;
3215: PetscFunctionReturn(PETSC_SUCCESS);
3216: }
3218: /*@
3219: PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element created with `PetscGaussLobattoLegendreElementAdvectionCreate()`
3221: Not Collective
3223: Input Parameters:
3224: + n - the number of GLL nodes
3225: . nodes - the GLL nodes, ignored
3226: . weights - the GLL weights, ignored
3227: - AA - advection obtained with `PetscGaussLobattoLegendreElementAdvectionCreate()`
3229: Level: beginner
3231: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementAdvectionCreate()`
3232: @*/
3233: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n, PetscReal nodes[], PetscReal weights[], PetscReal ***AA)
3234: {
3235: PetscFunctionBegin;
3236: PetscCall(PetscFree((*AA)[0]));
3237: PetscCall(PetscFree(*AA));
3238: *AA = NULL;
3239: PetscFunctionReturn(PETSC_SUCCESS);
3240: }
3242: /*@
3243: PetscGaussLobattoLegendreElementMassCreate - Build the elemental mass matrix for a single 1D Gauss-Lobatto-Legendre (GLL) spectral element
3245: Not Collective; No Fortran Support
3247: Input Parameters:
3248: + n - number of GLL nodes
3249: . nodes - the GLL quadrature nodes
3250: - weights - the GLL quadrature weights
3252: Output Parameter:
3253: . AA - newly allocated `n` x `n` mass matrix as `PetscReal **`
3255: Level: beginner
3257: Note:
3258: Free with `PetscGaussLobattoLegendreElementMassDestroy()`.
3260: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementMassDestroy()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionCreate()`
3261: @*/
3262: PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3263: {
3264: PetscReal **A;
3265: const PetscReal *gllweights = weights;
3266: const PetscInt glln = n;
3268: PetscFunctionBegin;
3269: PetscCall(PetscMalloc1(glln, &A));
3270: PetscCall(PetscMalloc1(glln * glln, &A[0]));
3271: for (PetscInt i = 1; i < glln; i++) A[i] = A[i - 1] + glln;
3272: if (glln == 1) A[0][0] = 0.;
3273: for (PetscInt i = 0; i < glln; i++) {
3274: for (PetscInt j = 0; j < glln; j++) {
3275: A[i][j] = 0.;
3276: if (j == i) A[i][j] = gllweights[i];
3277: }
3278: }
3279: *AA = A;
3280: PetscFunctionReturn(PETSC_SUCCESS);
3281: }
3283: /*@
3284: PetscGaussLobattoLegendreElementMassDestroy - Free a 1D GLL elemental mass matrix created with `PetscGaussLobattoLegendreElementMassCreate()`
3286: Not Collective; No Fortran Support
3288: Input Parameters:
3289: + n - number of GLL nodes (ignored)
3290: . nodes - the GLL quadrature nodes (ignored)
3291: . weights - the GLL quadrature weights (ignored)
3292: - AA - the mass matrix to free; `*AA` is set to `NULL` on return
3294: Level: beginner
3296: .seealso: `PetscGaussLobattoLegendreElementMassCreate()`, `PetscDTGaussLobattoLegendreQuadrature()`
3297: @*/
3298: PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3299: {
3300: PetscFunctionBegin;
3301: PetscCall(PetscFree((*AA)[0]));
3302: PetscCall(PetscFree(*AA));
3303: *AA = NULL;
3304: PetscFunctionReturn(PETSC_SUCCESS);
3305: }
3307: /*@
3308: PetscDTIndexToBary - convert an index into a barycentric coordinate.
3310: Input Parameters:
3311: + 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)
3312: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
3313: - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum)
3315: Output Parameter:
3316: . coord - will be filled with the barycentric coordinate, of length `n`
3318: Level: beginner
3320: Note:
3321: The indices map to barycentric coordinates in lexicographic order, where the first index is the
3322: least significant and the last index is the most significant.
3324: .seealso: `PetscDTBaryToIndex()`
3325: @*/
3326: PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[])
3327: {
3328: PetscInt c, d, s, total, subtotal, nexttotal;
3330: PetscFunctionBeginHot;
3331: PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
3332: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative");
3333: if (!len) {
3334: PetscCheck(!sum && !index, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
3335: PetscFunctionReturn(PETSC_SUCCESS);
3336: }
3337: for (c = 1, total = 1; c <= len; c++) {
3338: /* total is the number of ways to have a tuple of length c with sum */
3339: if (index < total) break;
3340: total = (total * (sum + c)) / c;
3341: }
3342: PetscCheck(c <= len, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range");
3343: for (d = c; d < len; d++) coord[d] = 0;
3344: for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) {
3345: /* subtotal is the number of ways to have a tuple of length c with sum s */
3346: /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */
3347: if (index + subtotal >= total) {
3348: coord[--c] = sum - s;
3349: index -= (total - subtotal);
3350: sum = s;
3351: total = nexttotal;
3352: subtotal = 1;
3353: nexttotal = 1;
3354: s = 0;
3355: } else {
3356: subtotal = (subtotal * (c + s)) / (s + 1);
3357: nexttotal = (nexttotal * (c - 1 + s)) / (s + 1);
3358: s++;
3359: }
3360: }
3361: PetscFunctionReturn(PETSC_SUCCESS);
3362: }
3364: /*@
3365: PetscDTBaryToIndex - convert a barycentric coordinate to an index
3367: Input Parameters:
3368: + 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)
3369: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
3370: - coord - a barycentric coordinate with the given length `len` and `sum`
3372: Output Parameter:
3373: . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum)
3375: Level: beginner
3377: Note:
3378: The indices map to barycentric coordinates in lexicographic order, where the first index is the
3379: least significant and the last index is the most significant.
3381: .seealso: `PetscDTIndexToBary`
3382: @*/
3383: PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index)
3384: {
3385: PetscInt c;
3386: PetscInt i;
3387: PetscInt total;
3389: PetscFunctionBeginHot;
3390: PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
3391: if (!len) {
3392: PetscCheck(!sum, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
3393: *index = 0;
3394: PetscFunctionReturn(PETSC_SUCCESS);
3395: }
3396: for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c;
3397: i = total - 1;
3398: c = len - 1;
3399: sum -= coord[c];
3400: while (sum > 0) {
3401: PetscInt subtotal;
3402: PetscInt s;
3404: for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s;
3405: i -= subtotal;
3406: sum -= coord[--c];
3407: }
3408: *index = i;
3409: PetscFunctionReturn(PETSC_SUCCESS);
3410: }
3412: /*@
3413: PetscQuadratureComputePermutations - Compute permutations of quadrature points corresponding to domain orientations
3415: Input Parameter:
3416: . quad - The `PetscQuadrature`
3418: Output Parameters:
3419: + Np - The number of domain orientations
3420: - perm - An array of `IS` permutations, one for ech orientation,
3422: Level: developer
3424: .seealso: `PetscQuadratureSetCellType()`, `PetscQuadrature`
3425: @*/
3426: PetscErrorCode PetscQuadratureComputePermutations(PetscQuadrature quad, PeOp PetscInt *Np, IS *perm[])
3427: {
3428: DMPolytopeType ct;
3429: const PetscReal *xq, *wq;
3430: PetscInt dim, qdim, d, Na, o, Nq, q, qp;
3432: PetscFunctionBegin;
3433: PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &xq, &wq));
3434: PetscCall(PetscQuadratureGetCellType(quad, &ct));
3435: dim = DMPolytopeTypeGetDim(ct);
3436: Na = DMPolytopeTypeGetNumArrangements(ct);
3437: PetscCall(PetscMalloc1(Na, perm));
3438: if (Np) *Np = Na;
3439: Na /= 2;
3440: for (o = -Na; o < Na; ++o) {
3441: DM refdm;
3442: PetscInt *idx;
3443: PetscReal xi0[3] = {-1., -1., -1.}, v0[3], J[9], detJ, txq[3];
3444: PetscBool flg;
3446: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
3447: PetscCall(DMPlexOrientPoint(refdm, 0, o));
3448: PetscCall(DMPlexComputeCellGeometryFEM(refdm, 0, NULL, v0, J, NULL, &detJ));
3449: PetscCall(PetscMalloc1(Nq, &idx));
3450: for (q = 0; q < Nq; ++q) {
3451: CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q * dim], txq);
3452: for (qp = 0; qp < Nq; ++qp) {
3453: PetscReal diff = 0.;
3455: for (d = 0; d < dim; ++d) diff += PetscAbsReal(txq[d] - xq[qp * dim + d]);
3456: if (diff < PETSC_SMALL) break;
3457: }
3458: PetscCheck(qp < Nq, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Transformed quad point %" PetscInt_FMT " does not match another quad point", q);
3459: idx[q] = qp;
3460: }
3461: PetscCall(DMDestroy(&refdm));
3462: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nq, idx, PETSC_OWN_POINTER, &(*perm)[o + Na]));
3463: PetscCall(ISGetInfo((*perm)[o + Na], IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg));
3464: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ordering for orientation %" PetscInt_FMT " was not a permutation", o);
3465: PetscCall(ISSetPermutation((*perm)[o + Na]));
3466: }
3467: if (!Na) (*perm)[0] = NULL;
3468: PetscFunctionReturn(PETSC_SUCCESS);
3469: }
3471: /*@
3472: PetscDTCreateQuadratureByCell - Create default quadrature for a given cell
3474: Not collective
3476: Input Parameters:
3477: + ct - The integration domain
3478: . qorder - The desired quadrature order
3479: - qtype - The type of simplex quadrature, or PETSCDTSIMPLEXQUAD_DEFAULT
3481: Output Parameters:
3482: + q - The cell quadrature
3483: - fq - The face quadrature
3485: Level: developer
3487: .seealso: `PetscDTCreateDefaultQuadrature()`, `PetscFECreateDefault()`, `PetscDTGaussTensorQuadrature()`, `PetscDTSimplexQuadrature()`, `PetscDTTensorQuadratureCreate()`
3488: @*/
3489: PetscErrorCode PetscDTCreateQuadratureByCell(DMPolytopeType ct, PetscInt qorder, PetscDTSimplexQuadratureType qtype, PetscQuadrature *q, PetscQuadrature *fq)
3490: {
3491: const PetscInt quadPointsPerEdge = PetscMax(qorder + 1, 1);
3492: const PetscInt dim = DMPolytopeTypeGetDim(ct);
3494: PetscFunctionBegin;
3495: switch (ct) {
3496: case DM_POLYTOPE_POINT: {
3497: PetscReal *x, *w;
3499: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
3500: PetscCall(PetscQuadratureSetCellType(*q, ct));
3501: PetscCall(PetscQuadratureSetOrder(*q, 0));
3502: PetscCall(PetscMalloc1(1, &x));
3503: PetscCall(PetscMalloc1(1, &w));
3504: x[0] = 0.;
3505: w[0] = 1.;
3506: PetscCall(PetscQuadratureSetData(*q, dim, 1, 1, x, w));
3507: *fq = NULL;
3508: } break;
3509: case DM_POLYTOPE_SEGMENT:
3510: case DM_POLYTOPE_POINT_PRISM_TENSOR:
3511: case DM_POLYTOPE_QUADRILATERAL:
3512: case DM_POLYTOPE_SEG_PRISM_TENSOR:
3513: case DM_POLYTOPE_HEXAHEDRON:
3514: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
3515: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q));
3516: PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq));
3517: break;
3518: case DM_POLYTOPE_TRIANGLE:
3519: case DM_POLYTOPE_TETRAHEDRON:
3520: PetscCall(PetscDTSimplexQuadrature(dim, 2 * qorder, qtype, q));
3521: PetscCall(PetscDTSimplexQuadrature(dim - 1, 2 * qorder, qtype, fq));
3522: break;
3523: case DM_POLYTOPE_TRI_PRISM:
3524: case DM_POLYTOPE_TRI_PRISM_TENSOR: {
3525: PetscQuadrature q1, q2;
3527: // TODO: this should be able to use symmetric rules, but doing so causes tests to fail
3528: PetscCall(PetscDTSimplexQuadrature(2, 2 * qorder, PETSCDTSIMPLEXQUAD_CONIC, &q1));
3529: PetscCall(PetscDTGaussTensorQuadrature(1, 1, quadPointsPerEdge, -1.0, 1.0, &q2));
3530: PetscCall(PetscDTTensorQuadratureCreate(q1, q2, q));
3531: PetscCall(PetscQuadratureDestroy(&q2));
3532: *fq = q1;
3533: /* TODO Need separate quadratures for each face */
3534: } break;
3535: default:
3536: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No quadrature for celltype %s", DMPolytopeTypes[PetscMin(ct, DM_POLYTOPE_UNKNOWN)]);
3537: }
3538: PetscFunctionReturn(PETSC_SUCCESS);
3539: }
3541: /*@
3542: PetscDTCreateDefaultQuadrature - Create default quadrature for a given cell
3544: Not collective
3546: Input Parameters:
3547: + ct - The integration domain
3548: - qorder - The desired quadrature order
3550: Output Parameters:
3551: + q - The cell quadrature
3552: - fq - The face quadrature
3554: Level: developer
3556: .seealso: `PetscDTCreateQuadratureByCell()`, `PetscFECreateDefault()`, `PetscDTGaussTensorQuadrature()`, `PetscDTSimplexQuadrature()`, `PetscDTTensorQuadratureCreate()`
3557: @*/
3558: PetscErrorCode PetscDTCreateDefaultQuadrature(DMPolytopeType ct, PetscInt qorder, PetscQuadrature *q, PetscQuadrature *fq)
3559: {
3560: PetscFunctionBegin;
3561: PetscCall(PetscDTCreateQuadratureByCell(ct, qorder, PETSCDTSIMPLEXQUAD_DEFAULT, q, fq));
3562: PetscFunctionReturn(PETSC_SUCCESS);
3563: }