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