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