Actual source code: dt.c

  1: /* Discretization tools */

  3: #include <petscdt.h>
  4: #include <petscblaslapack.h>
  5: #include <petsc/private/petscimpl.h>
  6: #include <petsc/private/dtimpl.h>
  7: #include <petsc/private/petscfeimpl.h>
  8: #include <petscviewer.h>
  9: #include <petscdmplex.h>
 10: #include <petscdmshell.h>

 12: #if defined(PETSC_HAVE_MPFR)
 13:   #include <mpfr.h>
 14: #endif

 16: const char *const        PetscDTNodeTypes_shifted[] = {"default", "gaussjacobi", "equispaced", "tanhsinh", "PETSCDTNODES_", NULL};
 17: const char *const *const PetscDTNodeTypes           = PetscDTNodeTypes_shifted + 1;

 19: const char *const        PetscDTSimplexQuadratureTypes_shifted[] = {"default", "conic", "minsym", "PETSCDTSIMPLEXQUAD_", NULL};
 20: const char *const *const PetscDTSimplexQuadratureTypes           = PetscDTSimplexQuadratureTypes_shifted + 1;

 22: static PetscBool GolubWelschCite       = PETSC_FALSE;
 23: const char       GolubWelschCitation[] = "@article{GolubWelsch1969,\n"
 24:                                          "  author  = {Golub and Welsch},\n"
 25:                                          "  title   = {Calculation of Quadrature Rules},\n"
 26:                                          "  journal = {Math. Comp.},\n"
 27:                                          "  volume  = {23},\n"
 28:                                          "  number  = {106},\n"
 29:                                          "  pages   = {221--230},\n"
 30:                                          "  year    = {1969}\n}\n";

 32: /* Numerical tests in src/dm/dt/tests/ex1.c show that when computing the nodes and weights of Gauss-Jacobi
 33:    quadrature rules:

 35:    - in double precision, Newton's method and Golub & Welsch both work for moderate degrees (< 100),
 36:    - in single precision, Newton's method starts producing incorrect roots around n = 15, but
 37:      the weights from Golub & Welsch become a problem before then: they produces errors
 38:      in computing the Jacobi-polynomial Gram matrix around n = 6.

 40:    So we default to Newton's method (required fewer dependencies) */
 41: PetscBool PetscDTGaussQuadratureNewton_Internal = PETSC_TRUE;

 43: PetscClassId PETSCQUADRATURE_CLASSID = 0;

 45: /*@
 46:   PetscQuadratureCreate - Create a `PetscQuadrature` object

 48:   Collective

 50:   Input Parameter:
 51: . comm - The communicator for the `PetscQuadrature` object

 53:   Output Parameter:
 54: . q  - The `PetscQuadrature` object

 56:   Level: beginner

 58: .seealso: `PetscQuadrature`, `Petscquadraturedestroy()`, `PetscQuadratureGetData()`
 59: @*/
 60: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
 61: {
 62:   PetscFunctionBegin;
 64:   PetscCall(DMInitializePackage());
 65:   PetscCall(PetscHeaderCreate(*q, PETSCQUADRATURE_CLASSID, "PetscQuadrature", "Quadrature", "DT", comm, PetscQuadratureDestroy, PetscQuadratureView));
 66:   (*q)->ct        = DM_POLYTOPE_UNKNOWN;
 67:   (*q)->dim       = -1;
 68:   (*q)->Nc        = 1;
 69:   (*q)->order     = -1;
 70:   (*q)->numPoints = 0;
 71:   (*q)->points    = NULL;
 72:   (*q)->weights   = NULL;
 73:   PetscFunctionReturn(PETSC_SUCCESS);
 74: }

 76: /*@
 77:   PetscQuadratureDuplicate - Create a deep copy of the `PetscQuadrature` object

 79:   Collective

 81:   Input Parameter:
 82: . q  - The `PetscQuadrature` object

 84:   Output Parameter:
 85: . r  - The new `PetscQuadrature` object

 87:   Level: beginner

 89: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureDestroy()`, `PetscQuadratureGetData()`
 90: @*/
 91: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
 92: {
 93:   DMPolytopeType   ct;
 94:   PetscInt         order, dim, Nc, Nq;
 95:   const PetscReal *points, *weights;
 96:   PetscReal       *p, *w;

 98:   PetscFunctionBegin;
100:   PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)q), r));
101:   PetscCall(PetscQuadratureGetCellType(q, &ct));
102:   PetscCall(PetscQuadratureSetCellType(*r, ct));
103:   PetscCall(PetscQuadratureGetOrder(q, &order));
104:   PetscCall(PetscQuadratureSetOrder(*r, order));
105:   PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights));
106:   PetscCall(PetscMalloc1(Nq * dim, &p));
107:   PetscCall(PetscMalloc1(Nq * Nc, &w));
108:   PetscCall(PetscArraycpy(p, points, Nq * dim));
109:   PetscCall(PetscArraycpy(w, weights, Nc * Nq));
110:   PetscCall(PetscQuadratureSetData(*r, dim, Nc, Nq, p, w));
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }

114: /*@
115:   PetscQuadratureDestroy - Destroys a `PetscQuadrature` object

117:   Collective

119:   Input Parameter:
120: . q  - The `PetscQuadrature` object

122:   Level: beginner

124: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
125: @*/
126: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
127: {
128:   PetscFunctionBegin;
129:   if (!*q) PetscFunctionReturn(PETSC_SUCCESS);
131:   if (--((PetscObject)(*q))->refct > 0) {
132:     *q = NULL;
133:     PetscFunctionReturn(PETSC_SUCCESS);
134:   }
135:   PetscCall(PetscFree((*q)->points));
136:   PetscCall(PetscFree((*q)->weights));
137:   PetscCall(PetscHeaderDestroy(q));
138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

141: /*@
142:   PetscQuadratureGetCellType - Return the cell type of the integration domain

144:   Not Collective

146:   Input Parameter:
147: . q - The `PetscQuadrature` object

149:   Output Parameter:
150: . ct - The cell type of the integration domain

152:   Level: intermediate

154: .seealso: `PetscQuadrature`, `PetscQuadratureSetCellType()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
155: @*/
156: PetscErrorCode PetscQuadratureGetCellType(PetscQuadrature q, DMPolytopeType *ct)
157: {
158:   PetscFunctionBegin;
161:   *ct = q->ct;
162:   PetscFunctionReturn(PETSC_SUCCESS);
163: }

165: /*@
166:   PetscQuadratureSetCellType - Set the cell type of the integration domain

168:   Not Collective

170:   Input Parameters:
171: + q - The `PetscQuadrature` object
172: - ct - The cell type of the integration domain

174:   Level: intermediate

176: .seealso: `PetscQuadrature`, `PetscQuadratureGetCellType()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
177: @*/
178: PetscErrorCode PetscQuadratureSetCellType(PetscQuadrature q, DMPolytopeType ct)
179: {
180:   PetscFunctionBegin;
182:   q->ct = ct;
183:   PetscFunctionReturn(PETSC_SUCCESS);
184: }

186: /*@
187:   PetscQuadratureGetOrder - Return the order of the method in the `PetscQuadrature`

189:   Not Collective

191:   Input Parameter:
192: . q - The `PetscQuadrature` object

194:   Output Parameter:
195: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated

197:   Level: intermediate

199: .seealso: `PetscQuadrature`, `PetscQuadratureSetOrder()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
200: @*/
201: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
202: {
203:   PetscFunctionBegin;
206:   *order = q->order;
207:   PetscFunctionReturn(PETSC_SUCCESS);
208: }

210: /*@
211:   PetscQuadratureSetOrder - Set the order of the method in the `PetscQuadrature`

213:   Not Collective

215:   Input Parameters:
216: + q - The `PetscQuadrature` object
217: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated

219:   Level: intermediate

221: .seealso: `PetscQuadrature`, `PetscQuadratureGetOrder()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
222: @*/
223: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
224: {
225:   PetscFunctionBegin;
227:   q->order = order;
228:   PetscFunctionReturn(PETSC_SUCCESS);
229: }

231: /*@
232:   PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated

234:   Not Collective

236:   Input Parameter:
237: . q - The `PetscQuadrature` object

239:   Output Parameter:
240: . Nc - The number of components

242:   Level: intermediate

244:   Note:
245:   We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.

247: .seealso: `PetscQuadrature`, `PetscQuadratureSetNumComponents()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
248: @*/
249: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
250: {
251:   PetscFunctionBegin;
254:   *Nc = q->Nc;
255:   PetscFunctionReturn(PETSC_SUCCESS);
256: }

258: /*@
259:   PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated

261:   Not Collective

263:   Input Parameters:
264: + q  - The `PetscQuadrature` object
265: - Nc - The number of components

267:   Level: intermediate

269:   Note:
270:   We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.

272: .seealso: `PetscQuadrature`, `PetscQuadratureGetNumComponents()`, `PetscQuadratureGetData()`, `PetscQuadratureSetData()`
273: @*/
274: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
275: {
276:   PetscFunctionBegin;
278:   q->Nc = Nc;
279:   PetscFunctionReturn(PETSC_SUCCESS);
280: }

282: /*@C
283:   PetscQuadratureGetData - Returns the data defining the `PetscQuadrature`

285:   Not Collective

287:   Input Parameter:
288: . q  - The `PetscQuadrature` object

290:   Output Parameters:
291: + dim - The spatial dimension
292: . Nc - The number of components
293: . npoints - The number of quadrature points
294: . points - The coordinates of each quadrature point
295: - weights - The weight of each quadrature point

297:   Level: intermediate

299:   Fortran Note:
300:   From Fortran you must call `PetscQuadratureRestoreData()` when you are done with the data

302: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureSetData()`
303: @*/
304: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
305: {
306:   PetscFunctionBegin;
308:   if (dim) {
310:     *dim = q->dim;
311:   }
312:   if (Nc) {
314:     *Nc = q->Nc;
315:   }
316:   if (npoints) {
318:     *npoints = q->numPoints;
319:   }
320:   if (points) {
322:     *points = q->points;
323:   }
324:   if (weights) {
326:     *weights = q->weights;
327:   }
328:   PetscFunctionReturn(PETSC_SUCCESS);
329: }

331: /*@
332:   PetscQuadratureEqual - determine whether two quadratures are equivalent

334:   Input Parameters:
335: + A - A `PetscQuadrature` object
336: - B - Another `PetscQuadrature` object

338:   Output Parameter:
339: . equal - `PETSC_TRUE` if the quadratures are the same

341:   Level: intermediate

343: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`
344: @*/
345: PetscErrorCode PetscQuadratureEqual(PetscQuadrature A, PetscQuadrature B, PetscBool *equal)
346: {
347:   PetscFunctionBegin;
351:   *equal = PETSC_FALSE;
352:   if (A->ct != B->ct || A->dim != B->dim || A->Nc != B->Nc || A->order != B->order || A->numPoints != B->numPoints) PetscFunctionReturn(PETSC_SUCCESS);
353:   for (PetscInt i = 0; i < A->numPoints * A->dim; i++) {
354:     if (PetscAbsReal(A->points[i] - B->points[i]) > PETSC_SMALL) PetscFunctionReturn(PETSC_SUCCESS);
355:   }
356:   if (!A->weights && !B->weights) {
357:     *equal = PETSC_TRUE;
358:     PetscFunctionReturn(PETSC_SUCCESS);
359:   }
360:   if (A->weights && B->weights) {
361:     for (PetscInt i = 0; i < A->numPoints; i++) {
362:       if (PetscAbsReal(A->weights[i] - B->weights[i]) > PETSC_SMALL) PetscFunctionReturn(PETSC_SUCCESS);
363:     }
364:     *equal = PETSC_TRUE;
365:   }
366:   PetscFunctionReturn(PETSC_SUCCESS);
367: }

369: static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[])
370: {
371:   PetscScalar *Js, *Jinvs;
372:   PetscInt     i, j, k;
373:   PetscBLASInt bm, bn, info;

375:   PetscFunctionBegin;
376:   if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
377:   PetscCall(PetscBLASIntCast(m, &bm));
378:   PetscCall(PetscBLASIntCast(n, &bn));
379: #if defined(PETSC_USE_COMPLEX)
380:   PetscCall(PetscMalloc2(m * n, &Js, m * n, &Jinvs));
381:   for (i = 0; i < m * n; i++) Js[i] = J[i];
382: #else
383:   Js    = (PetscReal *)J;
384:   Jinvs = Jinv;
385: #endif
386:   if (m == n) {
387:     PetscBLASInt *pivots;
388:     PetscScalar  *W;

390:     PetscCall(PetscMalloc2(m, &pivots, m, &W));

392:     PetscCall(PetscArraycpy(Jinvs, Js, m * m));
393:     PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info));
394:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscInt_FMT, (PetscInt)info);
395:     PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &info));
396:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscInt_FMT, (PetscInt)info);
397:     PetscCall(PetscFree2(pivots, W));
398:   } else if (m < n) {
399:     PetscScalar  *JJT;
400:     PetscBLASInt *pivots;
401:     PetscScalar  *W;

403:     PetscCall(PetscMalloc1(m * m, &JJT));
404:     PetscCall(PetscMalloc2(m, &pivots, m, &W));
405:     for (i = 0; i < m; i++) {
406:       for (j = 0; j < m; j++) {
407:         PetscScalar val = 0.;

409:         for (k = 0; k < n; k++) val += Js[i * n + k] * Js[j * n + k];
410:         JJT[i * m + j] = val;
411:       }
412:     }

414:     PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info));
415:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscInt_FMT, (PetscInt)info);
416:     PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info));
417:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscInt_FMT, (PetscInt)info);
418:     for (i = 0; i < n; i++) {
419:       for (j = 0; j < m; j++) {
420:         PetscScalar val = 0.;

422:         for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j];
423:         Jinvs[i * m + j] = val;
424:       }
425:     }
426:     PetscCall(PetscFree2(pivots, W));
427:     PetscCall(PetscFree(JJT));
428:   } else {
429:     PetscScalar  *JTJ;
430:     PetscBLASInt *pivots;
431:     PetscScalar  *W;

433:     PetscCall(PetscMalloc1(n * n, &JTJ));
434:     PetscCall(PetscMalloc2(n, &pivots, n, &W));
435:     for (i = 0; i < n; i++) {
436:       for (j = 0; j < n; j++) {
437:         PetscScalar val = 0.;

439:         for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j];
440:         JTJ[i * n + j] = val;
441:       }
442:     }

444:     PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bn, pivots, &info));
445:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscInt_FMT, (PetscInt)info);
446:     PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info));
447:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscInt_FMT, (PetscInt)info);
448:     for (i = 0; i < n; i++) {
449:       for (j = 0; j < m; j++) {
450:         PetscScalar val = 0.;

452:         for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k];
453:         Jinvs[i * m + j] = val;
454:       }
455:     }
456:     PetscCall(PetscFree2(pivots, W));
457:     PetscCall(PetscFree(JTJ));
458:   }
459: #if defined(PETSC_USE_COMPLEX)
460:   for (i = 0; i < m * n; i++) Jinv[i] = PetscRealPart(Jinvs[i]);
461:   PetscCall(PetscFree2(Js, Jinvs));
462: #endif
463:   PetscFunctionReturn(PETSC_SUCCESS);
464: }

466: /*@
467:    PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation.

469:    Collective

471:    Input Parameters:
472: +  q - the quadrature functional
473: .  imageDim - the dimension of the image of the transformation
474: .  origin - a point in the original space
475: .  originImage - the image of the origin under the transformation
476: .  J - the Jacobian of the image: an [imageDim x dim] matrix in row major order
477: -  formDegree - transform the quadrature weights as k-forms of this form degree (if the number of components is a multiple of (dim choose formDegree), it is assumed that they represent multiple k-forms) [see `PetscDTAltVPullback()` for interpretation of formDegree]

479:    Output Parameter:
480: .  Jinvstarq - a quadrature rule where each point is the image of a point in the original quadrature rule, and where the k-form weights have been pulled-back by the pseudoinverse of `J` to the k-form weights in the image space.

482:    Level: intermediate

484:    Note:
485:    The new quadrature rule will have a different number of components if spaces have different dimensions.  For example, pushing a 2-form forward from a two dimensional space to a three dimensional space changes the number of components from 1 to 3.

487: .seealso: `PetscQuadrature`, `PetscDTAltVPullback()`, `PetscDTAltVPullbackMatrix()`
488: @*/
489: PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq)
490: {
491:   PetscInt         dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c;
492:   const PetscReal *points;
493:   const PetscReal *weights;
494:   PetscReal       *imagePoints, *imageWeights;
495:   PetscReal       *Jinv;
496:   PetscReal       *Jinvstar;

498:   PetscFunctionBegin;
500:   PetscCheck(imageDim >= PetscAbsInt(formDegree), PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Cannot represent a %" PetscInt_FMT "-form in %" PetscInt_FMT " dimensions", PetscAbsInt(formDegree), imageDim);
501:   PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights));
502:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize));
503:   PetscCheck(Nc % formSize == 0, PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Number of components %" PetscInt_FMT " is not a multiple of formSize %" PetscInt_FMT, Nc, formSize);
504:   Ncopies = Nc / formSize;
505:   PetscCall(PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize));
506:   imageNc = Ncopies * imageFormSize;
507:   PetscCall(PetscMalloc1(Npoints * imageDim, &imagePoints));
508:   PetscCall(PetscMalloc1(Npoints * imageNc, &imageWeights));
509:   PetscCall(PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar));
510:   PetscCall(PetscDTJacobianInverse_Internal(imageDim, dim, J, Jinv));
511:   PetscCall(PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar));
512:   for (pt = 0; pt < Npoints; pt++) {
513:     const PetscReal *point      = &points[pt * dim];
514:     PetscReal       *imagePoint = &imagePoints[pt * imageDim];

516:     for (i = 0; i < imageDim; i++) {
517:       PetscReal val = originImage[i];

519:       for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]);
520:       imagePoint[i] = val;
521:     }
522:     for (c = 0; c < Ncopies; c++) {
523:       const PetscReal *form      = &weights[pt * Nc + c * formSize];
524:       PetscReal       *imageForm = &imageWeights[pt * imageNc + c * imageFormSize];

526:       for (i = 0; i < imageFormSize; i++) {
527:         PetscReal val = 0.;

529:         for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j];
530:         imageForm[i] = val;
531:       }
532:     }
533:   }
534:   PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq));
535:   PetscCall(PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights));
536:   PetscCall(PetscFree2(Jinv, Jinvstar));
537:   PetscFunctionReturn(PETSC_SUCCESS);
538: }

540: /*@C
541:   PetscQuadratureSetData - Sets the data defining the quadrature

543:   Not Collective

545:   Input Parameters:
546: + q  - The `PetscQuadrature` object
547: . dim - The spatial dimension
548: . Nc - The number of components
549: . npoints - The number of quadrature points
550: . points - The coordinates of each quadrature point
551: - weights - The weight of each quadrature point

553:   Level: intermediate

555:   Note:
556:   This routine owns the references to points and weights, so they must be allocated using `PetscMalloc()` and the user should not free them.

558: .seealso: `PetscQuadrature`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
559: @*/
560: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
561: {
562:   PetscFunctionBegin;
564:   if (dim >= 0) q->dim = dim;
565:   if (Nc >= 0) q->Nc = Nc;
566:   if (npoints >= 0) q->numPoints = npoints;
567:   if (points) {
569:     q->points = points;
570:   }
571:   if (weights) {
573:     q->weights = weights;
574:   }
575:   PetscFunctionReturn(PETSC_SUCCESS);
576: }

578: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
579: {
580:   PetscInt          q, d, c;
581:   PetscViewerFormat format;

583:   PetscFunctionBegin;
584:   if (quad->Nc > 1)
585:     PetscCall(PetscViewerASCIIPrintf(v, "Quadrature on a %s of order %" PetscInt_FMT " on %" PetscInt_FMT " points (dim %" PetscInt_FMT ") with %" PetscInt_FMT " components\n", DMPolytopeTypes[quad->ct], quad->order, quad->numPoints, quad->dim, quad->Nc));
586:   else PetscCall(PetscViewerASCIIPrintf(v, "Quadrature on a %s of order %" PetscInt_FMT " on %" PetscInt_FMT " points (dim %" PetscInt_FMT ")\n", DMPolytopeTypes[quad->ct], quad->order, quad->numPoints, quad->dim));
587:   PetscCall(PetscViewerGetFormat(v, &format));
588:   if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
589:   for (q = 0; q < quad->numPoints; ++q) {
590:     PetscCall(PetscViewerASCIIPrintf(v, "p%" PetscInt_FMT " (", q));
591:     PetscCall(PetscViewerASCIIUseTabs(v, PETSC_FALSE));
592:     for (d = 0; d < quad->dim; ++d) {
593:       if (d) PetscCall(PetscViewerASCIIPrintf(v, ", "));
594:       PetscCall(PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q * quad->dim + d]));
595:     }
596:     PetscCall(PetscViewerASCIIPrintf(v, ") "));
597:     if (quad->Nc > 1) PetscCall(PetscViewerASCIIPrintf(v, "w%" PetscInt_FMT " (", q));
598:     for (c = 0; c < quad->Nc; ++c) {
599:       if (c) PetscCall(PetscViewerASCIIPrintf(v, ", "));
600:       PetscCall(PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q * quad->Nc + c]));
601:     }
602:     if (quad->Nc > 1) PetscCall(PetscViewerASCIIPrintf(v, ")"));
603:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
604:     PetscCall(PetscViewerASCIIUseTabs(v, PETSC_TRUE));
605:   }
606:   PetscFunctionReturn(PETSC_SUCCESS);
607: }

609: /*@C
610:   PetscQuadratureView - View a `PetscQuadrature` object

612:   Collective

614:   Input Parameters:
615: + quad  - The `PetscQuadrature` object
616: - viewer - The `PetscViewer` object

618:   Level: beginner

620: .seealso: `PetscQuadrature`, `PetscViewer`, `PetscQuadratureCreate()`, `PetscQuadratureGetData()`
621: @*/
622: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
623: {
624:   PetscBool iascii;

626:   PetscFunctionBegin;
629:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)quad), &viewer));
630:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
631:   PetscCall(PetscViewerASCIIPushTab(viewer));
632:   if (iascii) PetscCall(PetscQuadratureView_Ascii(quad, viewer));
633:   PetscCall(PetscViewerASCIIPopTab(viewer));
634:   PetscFunctionReturn(PETSC_SUCCESS);
635: }

637: /*@C
638:   PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement

640:   Not Collective; No Fortran Support

642:   Input Parameters:
643: + q - The original `PetscQuadrature`
644: . numSubelements - The number of subelements the original element is divided into
645: . v0 - An array of the initial points for each subelement
646: - jac - An array of the Jacobian mappings from the reference to each subelement

648:   Output Parameter:
649: . dim - The dimension

651:   Level: intermediate

653:   Note:
654:   Together v0 and jac define an affine mapping from the original reference element to each subelement

656: .seealso: `PetscQuadrature`, `PetscFECreate()`, `PetscSpaceGetDimension()`, `PetscDualSpaceGetDimension()`
657: @*/
658: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
659: {
660:   DMPolytopeType   ct;
661:   const PetscReal *points, *weights;
662:   PetscReal       *pointsRef, *weightsRef;
663:   PetscInt         dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;

665:   PetscFunctionBegin;
670:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, qref));
671:   PetscCall(PetscQuadratureGetCellType(q, &ct));
672:   PetscCall(PetscQuadratureGetOrder(q, &order));
673:   PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
674:   npointsRef = npoints * numSubelements;
675:   PetscCall(PetscMalloc1(npointsRef * dim, &pointsRef));
676:   PetscCall(PetscMalloc1(npointsRef * Nc, &weightsRef));
677:   for (c = 0; c < numSubelements; ++c) {
678:     for (p = 0; p < npoints; ++p) {
679:       for (d = 0; d < dim; ++d) {
680:         pointsRef[(c * npoints + p) * dim + d] = v0[c * dim + d];
681:         for (e = 0; e < dim; ++e) pointsRef[(c * npoints + p) * dim + d] += jac[(c * dim + d) * dim + e] * (points[p * dim + e] + 1.0);
682:       }
683:       /* Could also use detJ here */
684:       for (cp = 0; cp < Nc; ++cp) weightsRef[(c * npoints + p) * Nc + cp] = weights[p * Nc + cp] / numSubelements;
685:     }
686:   }
687:   PetscCall(PetscQuadratureSetCellType(*qref, ct));
688:   PetscCall(PetscQuadratureSetOrder(*qref, order));
689:   PetscCall(PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef));
690:   PetscFunctionReturn(PETSC_SUCCESS);
691: }

693: /* Compute the coefficients for the Jacobi polynomial recurrence,
694:  *
695:  * J^{a,b}_n(x) = (cnm1 + cnm1x * x) * J^{a,b}_{n-1}(x) - cnm2 * J^{a,b}_{n-2}(x).
696:  */
697: #define PetscDTJacobiRecurrence_Internal(n, a, b, cnm1, cnm1x, cnm2) \
698:   do { \
699:     PetscReal _a = (a); \
700:     PetscReal _b = (b); \
701:     PetscReal _n = (n); \
702:     if (n == 1) { \
703:       (cnm1)  = (_a - _b) * 0.5; \
704:       (cnm1x) = (_a + _b + 2.) * 0.5; \
705:       (cnm2)  = 0.; \
706:     } else { \
707:       PetscReal _2n  = _n + _n; \
708:       PetscReal _d   = (_2n * (_n + _a + _b) * (_2n + _a + _b - 2)); \
709:       PetscReal _n1  = (_2n + _a + _b - 1.) * (_a * _a - _b * _b); \
710:       PetscReal _n1x = (_2n + _a + _b - 1.) * (_2n + _a + _b) * (_2n + _a + _b - 2); \
711:       PetscReal _n2  = 2. * ((_n + _a - 1.) * (_n + _b - 1.) * (_2n + _a + _b)); \
712:       (cnm1)         = _n1 / _d; \
713:       (cnm1x)        = _n1x / _d; \
714:       (cnm2)         = _n2 / _d; \
715:     } \
716:   } while (0)

718: /*@
719:   PetscDTJacobiNorm - Compute the weighted L2 norm of a Jacobi polynomial.

721:   $\| P^{\alpha,\beta}_n \|_{\alpha,\beta}^2 = \int_{-1}^1 (1 + x)^{\alpha} (1 - x)^{\beta} P^{\alpha,\beta}_n (x)^2 dx.$

723:   Input Parameters:
724: - alpha - the left exponent > -1
725: . beta - the right exponent > -1
726: + n - the polynomial degree

728:   Output Parameter:
729: . norm - the weighted L2 norm

731:   Level: beginner

733: .seealso: `PetscQuadrature`, `PetscDTJacobiEval()`
734: @*/
735: PetscErrorCode PetscDTJacobiNorm(PetscReal alpha, PetscReal beta, PetscInt n, PetscReal *norm)
736: {
737:   PetscReal twoab1;
738:   PetscReal gr;

740:   PetscFunctionBegin;
741:   PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent alpha %g <= -1. invalid", (double)alpha);
742:   PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent beta %g <= -1. invalid", (double)beta);
743:   PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "n %" PetscInt_FMT " < 0 invalid", n);
744:   twoab1 = PetscPowReal(2., alpha + beta + 1.);
745: #if defined(PETSC_HAVE_LGAMMA)
746:   if (!n) {
747:     gr = PetscExpReal(PetscLGamma(alpha + 1.) + PetscLGamma(beta + 1.) - PetscLGamma(alpha + beta + 2.));
748:   } else {
749:     gr = PetscExpReal(PetscLGamma(n + alpha + 1.) + PetscLGamma(n + beta + 1.) - (PetscLGamma(n + 1.) + PetscLGamma(n + alpha + beta + 1.))) / (n + n + alpha + beta + 1.);
750:   }
751: #else
752:   {
753:     PetscInt alphai = (PetscInt)alpha;
754:     PetscInt betai  = (PetscInt)beta;
755:     PetscInt i;

757:     gr = n ? (1. / (n + n + alpha + beta + 1.)) : 1.;
758:     if ((PetscReal)alphai == alpha) {
759:       if (!n) {
760:         for (i = 0; i < alphai; i++) gr *= (i + 1.) / (beta + i + 1.);
761:         gr /= (alpha + beta + 1.);
762:       } else {
763:         for (i = 0; i < alphai; i++) gr *= (n + i + 1.) / (n + beta + i + 1.);
764:       }
765:     } else if ((PetscReal)betai == beta) {
766:       if (!n) {
767:         for (i = 0; i < betai; i++) gr *= (i + 1.) / (alpha + i + 2.);
768:         gr /= (alpha + beta + 1.);
769:       } else {
770:         for (i = 0; i < betai; i++) gr *= (n + i + 1.) / (n + alpha + i + 1.);
771:       }
772:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
773:   }
774: #endif
775:   *norm = PetscSqrtReal(twoab1 * gr);
776:   PetscFunctionReturn(PETSC_SUCCESS);
777: }

779: static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p)
780: {
781:   PetscReal ak, bk;
782:   PetscReal abk1;
783:   PetscInt  i, l, maxdegree;

785:   PetscFunctionBegin;
786:   maxdegree = degrees[ndegree - 1] - k;
787:   ak        = a + k;
788:   bk        = b + k;
789:   abk1      = a + b + k + 1.;
790:   if (maxdegree < 0) {
791:     for (i = 0; i < npoints; i++)
792:       for (l = 0; l < ndegree; l++) p[i * ndegree + l] = 0.;
793:     PetscFunctionReturn(PETSC_SUCCESS);
794:   }
795:   for (i = 0; i < npoints; i++) {
796:     PetscReal pm1, pm2, x;
797:     PetscReal cnm1, cnm1x, cnm2;
798:     PetscInt  j, m;

800:     x   = points[i];
801:     pm2 = 1.;
802:     PetscDTJacobiRecurrence_Internal(1, ak, bk, cnm1, cnm1x, cnm2);
803:     pm1 = (cnm1 + cnm1x * x);
804:     l   = 0;
805:     while (l < ndegree && degrees[l] - k < 0) p[l++] = 0.;
806:     while (l < ndegree && degrees[l] - k == 0) {
807:       p[l] = pm2;
808:       for (m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5;
809:       l++;
810:     }
811:     while (l < ndegree && degrees[l] - k == 1) {
812:       p[l] = pm1;
813:       for (m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5;
814:       l++;
815:     }
816:     for (j = 2; j <= maxdegree; j++) {
817:       PetscReal pp;

819:       PetscDTJacobiRecurrence_Internal(j, ak, bk, cnm1, cnm1x, cnm2);
820:       pp  = (cnm1 + cnm1x * x) * pm1 - cnm2 * pm2;
821:       pm2 = pm1;
822:       pm1 = pp;
823:       while (l < ndegree && degrees[l] - k == j) {
824:         p[l] = pp;
825:         for (m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5;
826:         l++;
827:       }
828:     }
829:     p += ndegree;
830:   }
831:   PetscFunctionReturn(PETSC_SUCCESS);
832: }

834: /*@
835:   PetscDTJacobiEvalJet - Evaluate the jet (function and derivatives) of the Jacobi polynomials basis up to a given degree.
836:   The Jacobi polynomials with indices $\alpha$ and $\beta$ are orthogonal with respect to the weighted inner product
837:   $\langle f, g \rangle = \int_{-1}^1 (1+x)^{\alpha} (1-x)^{\beta} f(x) g(x) dx$.

839:   Input Parameters:
840: + alpha - the left exponent of the weight
841: . beta - the right exponetn of the weight
842: . npoints - the number of points to evaluate the polynomials at
843: . points - [npoints] array of point coordinates
844: . degree - the maximm degree polynomial space to evaluate, (degree + 1) will be evaluated total.
845: - k - the maximum derivative to evaluate in the jet, (k + 1) will be evaluated total.

847:   Output Parameter:
848: . p - an array containing the evaluations of the Jacobi polynomials's jets on the points.  the size is (degree + 1) x
849:   (k + 1) x npoints, which also describes the order of the dimensions of this three-dimensional array: the first
850:   (slowest varying) dimension is polynomial degree; the second dimension is derivative order; the third (fastest
851:   varying) dimension is the index of the evaluation point.

853:   Level: advanced

855: .seealso: `PetscDTJacobiEval()`, `PetscDTPKDEvalJet()`
856: @*/
857: PetscErrorCode PetscDTJacobiEvalJet(PetscReal alpha, PetscReal beta, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
858: {
859:   PetscInt   i, j, l;
860:   PetscInt  *degrees;
861:   PetscReal *psingle;

863:   PetscFunctionBegin;
864:   if (degree == 0) {
865:     PetscInt zero = 0;

867:     for (i = 0; i <= k; i++) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, 1, &zero, &p[i * npoints]));
868:     PetscFunctionReturn(PETSC_SUCCESS);
869:   }
870:   PetscCall(PetscMalloc1(degree + 1, &degrees));
871:   PetscCall(PetscMalloc1((degree + 1) * npoints, &psingle));
872:   for (i = 0; i <= degree; i++) degrees[i] = i;
873:   for (i = 0; i <= k; i++) {
874:     PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, degree + 1, degrees, psingle));
875:     for (j = 0; j <= degree; j++) {
876:       for (l = 0; l < npoints; l++) p[(j * (k + 1) + i) * npoints + l] = psingle[l * (degree + 1) + j];
877:     }
878:   }
879:   PetscCall(PetscFree(psingle));
880:   PetscCall(PetscFree(degrees));
881:   PetscFunctionReturn(PETSC_SUCCESS);
882: }

884: /*@
885:    PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$ at a set of points
886:                        at points

888:    Not Collective

890:    Input Parameters:
891: +  npoints - number of spatial points to evaluate at
892: .  alpha - the left exponent > -1
893: .  beta - the right exponent > -1
894: .  points - array of locations to evaluate at
895: .  ndegree - number of basis degrees to evaluate
896: -  degrees - sorted array of degrees to evaluate

898:    Output Parameters:
899: +  B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
900: .  D - row-oriented derivative evaluation matrix (or NULL)
901: -  D2 - row-oriented second derivative evaluation matrix (or NULL)

903:    Level: intermediate

905: .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()`
906: @*/
907: PetscErrorCode PetscDTJacobiEval(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *B, PetscReal *D, PetscReal *D2)
908: {
909:   PetscFunctionBegin;
910:   PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1.");
911:   PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1.");
912:   if (!npoints || !ndegree) PetscFunctionReturn(PETSC_SUCCESS);
913:   if (B) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B));
914:   if (D) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D));
915:   if (D2) PetscCall(PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2));
916:   PetscFunctionReturn(PETSC_SUCCESS);
917: }

919: /*@
920:    PetscDTLegendreEval - evaluate Legendre polynomials at points

922:    Not Collective

924:    Input Parameters:
925: +  npoints - number of spatial points to evaluate at
926: .  points - array of locations to evaluate at
927: .  ndegree - number of basis degrees to evaluate
928: -  degrees - sorted array of degrees to evaluate

930:    Output Parameters:
931: +  B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
932: .  D - row-oriented derivative evaluation matrix (or NULL)
933: -  D2 - row-oriented second derivative evaluation matrix (or NULL)

935:    Level: intermediate

937: .seealso: `PetscDTGaussQuadrature()`
938: @*/
939: PetscErrorCode PetscDTLegendreEval(PetscInt npoints, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *B, PetscReal *D, PetscReal *D2)
940: {
941:   PetscFunctionBegin;
942:   PetscCall(PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2));
943:   PetscFunctionReturn(PETSC_SUCCESS);
944: }

946: /*@
947:   PetscDTIndexToGradedOrder - convert an index into a tuple of monomial degrees in a graded order (that is, if the degree sum of tuple x is less than the degree sum of tuple y, then the index of x is smaller than the index of y)

949:   Input Parameters:
950: + len - the desired length of the degree tuple
951: - index - the index to convert: should be >= 0

953:   Output Parameter:
954: . degtup - will be filled with a tuple of degrees

956:   Level: beginner

958:   Note:
959:   For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
960:   acts as a tiebreaker.  For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
961:   last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).

963: .seealso: `PetscDTGradedOrderToIndex()`
964: @*/
965: PetscErrorCode PetscDTIndexToGradedOrder(PetscInt len, PetscInt index, PetscInt degtup[])
966: {
967:   PetscInt i, total;
968:   PetscInt sum;

970:   PetscFunctionBeginHot;
971:   PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
972:   PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative");
973:   total = 1;
974:   sum   = 0;
975:   while (index >= total) {
976:     index -= total;
977:     total = (total * (len + sum)) / (sum + 1);
978:     sum++;
979:   }
980:   for (i = 0; i < len; i++) {
981:     PetscInt c;

983:     degtup[i] = sum;
984:     for (c = 0, total = 1; c < sum; c++) {
985:       /* going into the loop, total is the number of way to have a tuple of sum exactly c with length len - 1 - i */
986:       if (index < total) break;
987:       index -= total;
988:       total = (total * (len - 1 - i + c)) / (c + 1);
989:       degtup[i]--;
990:     }
991:     sum -= degtup[i];
992:   }
993:   PetscFunctionReturn(PETSC_SUCCESS);
994: }

996: /*@
997:   PetscDTGradedOrderToIndex - convert a tuple into an index in a graded order, the inverse of `PetscDTIndexToGradedOrder()`.

999:   Input Parameters:
1000: + len - the length of the degree tuple
1001: - degtup - tuple with this length

1003:   Output Parameter:
1004: . index - index in graded order: >= 0

1006:   Level: Beginner

1008:   Note:
1009:   For two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
1010:   acts as a tiebreaker.  For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
1011:   last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).

1013: .seealso: `PetscDTIndexToGradedOrder()`
1014: @*/
1015: PetscErrorCode PetscDTGradedOrderToIndex(PetscInt len, const PetscInt degtup[], PetscInt *index)
1016: {
1017:   PetscInt i, idx, sum, total;

1019:   PetscFunctionBeginHot;
1020:   PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
1021:   for (i = 0, sum = 0; i < len; i++) sum += degtup[i];
1022:   idx   = 0;
1023:   total = 1;
1024:   for (i = 0; i < sum; i++) {
1025:     idx += total;
1026:     total = (total * (len + i)) / (i + 1);
1027:   }
1028:   for (i = 0; i < len - 1; i++) {
1029:     PetscInt c;

1031:     total = 1;
1032:     sum -= degtup[i];
1033:     for (c = 0; c < sum; c++) {
1034:       idx += total;
1035:       total = (total * (len - 1 - i + c)) / (c + 1);
1036:     }
1037:   }
1038:   *index = idx;
1039:   PetscFunctionReturn(PETSC_SUCCESS);
1040: }

1042: static PetscBool PKDCite       = PETSC_FALSE;
1043: const char       PKDCitation[] = "@article{Kirby2010,\n"
1044:                                  "  title={Singularity-free evaluation of collapsed-coordinate orthogonal polynomials},\n"
1045:                                  "  author={Kirby, Robert C},\n"
1046:                                  "  journal={ACM Transactions on Mathematical Software (TOMS)},\n"
1047:                                  "  volume={37},\n"
1048:                                  "  number={1},\n"
1049:                                  "  pages={1--16},\n"
1050:                                  "  year={2010},\n"
1051:                                  "  publisher={ACM New York, NY, USA}\n}\n";

1053: /*@
1054:   PetscDTPKDEvalJet - Evaluate the jet (function and derivatives) of the Proriol-Koornwinder-Dubiner (PKD) basis for
1055:   the space of polynomials up to a given degree.  The PKD basis is L2-orthonormal on the biunit simplex (which is used
1056:   as the reference element for finite elements in PETSc), which makes it a stable basis to use for evaluating
1057:   polynomials in that domain.

1059:   Input Parameters:
1060: + dim - the number of variables in the multivariate polynomials
1061: . npoints - the number of points to evaluate the polynomials at
1062: . points - [npoints x dim] array of point coordinates
1063: . degree - the degree (sum of degrees on the variables in a monomial) of the polynomial space to evaluate.  There are ((dim + degree) choose dim) polynomials in this space.
1064: - k - the maximum order partial derivative to evaluate in the jet.  There are (dim + k choose dim) partial derivatives
1065:   in the jet.  Choosing k = 0 means to evaluate just the function and no derivatives

1067:   Output Parameter:
1068: . p - an array containing the evaluations of the PKD polynomials' jets on the points.  The size is ((dim + degree)
1069:   choose dim) x ((dim + k) choose dim) x npoints, which also describes the order of the dimensions of this
1070:   three-dimensional array: the first (slowest varying) dimension is basis function index; the second dimension is jet
1071:   index; the third (fastest varying) dimension is the index of the evaluation point.

1073:   Level: advanced

1075:   Notes:
1076:   The ordering of the basis functions, and the ordering of the derivatives in the jet, both follow the graded
1077:   ordering of `PetscDTIndexToGradedOrder()` and `PetscDTGradedOrderToIndex()`.  For example, in 3D, the polynomial with
1078:   leading monomial x^2,y^0,z^1, which has degree tuple (2,0,1), which by `PetscDTGradedOrderToIndex()` has index 12 (it is the 13th basis function in the space);
1079:   the partial derivative $\partial_x \partial_z$ has order tuple (1,0,1), appears at index 6 in the jet (it is the 7th partial derivative in the jet).

1081:   The implementation uses Kirby's singularity-free evaluation algorithm, https://doi.org/10.1145/1644001.1644006.

1083: .seealso: `PetscDTGradedOrderToIndex()`, `PetscDTIndexToGradedOrder()`, `PetscDTJacobiEvalJet()`
1084: @*/
1085: PetscErrorCode PetscDTPKDEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
1086: {
1087:   PetscInt   degidx, kidx, d, pt;
1088:   PetscInt   Nk, Ndeg;
1089:   PetscInt  *ktup, *degtup;
1090:   PetscReal *scales, initscale, scaleexp;

1092:   PetscFunctionBegin;
1093:   PetscCall(PetscCitationsRegister(PKDCitation, &PKDCite));
1094:   PetscCall(PetscDTBinomialInt(dim + k, k, &Nk));
1095:   PetscCall(PetscDTBinomialInt(degree + dim, degree, &Ndeg));
1096:   PetscCall(PetscMalloc2(dim, &degtup, dim, &ktup));
1097:   PetscCall(PetscMalloc1(Ndeg, &scales));
1098:   initscale = 1.;
1099:   if (dim > 1) {
1100:     PetscCall(PetscDTBinomial(dim, 2, &scaleexp));
1101:     initscale = PetscPowReal(2., scaleexp * 0.5);
1102:   }
1103:   for (degidx = 0; degidx < Ndeg; degidx++) {
1104:     PetscInt  e, i;
1105:     PetscInt  m1idx = -1, m2idx = -1;
1106:     PetscInt  n;
1107:     PetscInt  degsum;
1108:     PetscReal alpha;
1109:     PetscReal cnm1, cnm1x, cnm2;
1110:     PetscReal norm;

1112:     PetscCall(PetscDTIndexToGradedOrder(dim, degidx, degtup));
1113:     for (d = dim - 1; d >= 0; d--)
1114:       if (degtup[d]) break;
1115:     if (d < 0) { /* constant is 1 everywhere, all derivatives are zero */
1116:       scales[degidx] = initscale;
1117:       for (e = 0; e < dim; e++) {
1118:         PetscCall(PetscDTJacobiNorm(e, 0., 0, &norm));
1119:         scales[degidx] /= norm;
1120:       }
1121:       for (i = 0; i < npoints; i++) p[degidx * Nk * npoints + i] = 1.;
1122:       for (i = 0; i < (Nk - 1) * npoints; i++) p[(degidx * Nk + 1) * npoints + i] = 0.;
1123:       continue;
1124:     }
1125:     n = degtup[d];
1126:     degtup[d]--;
1127:     PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m1idx));
1128:     if (degtup[d] > 0) {
1129:       degtup[d]--;
1130:       PetscCall(PetscDTGradedOrderToIndex(dim, degtup, &m2idx));
1131:       degtup[d]++;
1132:     }
1133:     degtup[d]++;
1134:     for (e = 0, degsum = 0; e < d; e++) degsum += degtup[e];
1135:     alpha = 2 * degsum + d;
1136:     PetscDTJacobiRecurrence_Internal(n, alpha, 0., cnm1, cnm1x, cnm2);

1138:     scales[degidx] = initscale;
1139:     for (e = 0, degsum = 0; e < dim; e++) {
1140:       PetscInt  f;
1141:       PetscReal ealpha;
1142:       PetscReal enorm;

1144:       ealpha = 2 * degsum + e;
1145:       for (f = 0; f < degsum; f++) scales[degidx] *= 2.;
1146:       PetscCall(PetscDTJacobiNorm(ealpha, 0., degtup[e], &enorm));
1147:       scales[degidx] /= enorm;
1148:       degsum += degtup[e];
1149:     }

1151:     for (pt = 0; pt < npoints; pt++) {
1152:       /* compute the multipliers */
1153:       PetscReal thetanm1, thetanm1x, thetanm2;

1155:       thetanm1x = dim - (d + 1) + 2. * points[pt * dim + d];
1156:       for (e = d + 1; e < dim; e++) thetanm1x += points[pt * dim + e];
1157:       thetanm1x *= 0.5;
1158:       thetanm1 = (2. - (dim - (d + 1)));
1159:       for (e = d + 1; e < dim; e++) thetanm1 -= points[pt * dim + e];
1160:       thetanm1 *= 0.5;
1161:       thetanm2 = thetanm1 * thetanm1;

1163:       for (kidx = 0; kidx < Nk; kidx++) {
1164:         PetscInt f;

1166:         PetscCall(PetscDTIndexToGradedOrder(dim, kidx, ktup));
1167:         /* first sum in the same derivative terms */
1168:         p[(degidx * Nk + kidx) * npoints + pt] = (cnm1 * thetanm1 + cnm1x * thetanm1x) * p[(m1idx * Nk + kidx) * npoints + pt];
1169:         if (m2idx >= 0) p[(degidx * Nk + kidx) * npoints + pt] -= cnm2 * thetanm2 * p[(m2idx * Nk + kidx) * npoints + pt];

1171:         for (f = d; f < dim; f++) {
1172:           PetscInt km1idx, mplty = ktup[f];

1174:           if (!mplty) continue;
1175:           ktup[f]--;
1176:           PetscCall(PetscDTGradedOrderToIndex(dim, ktup, &km1idx));

1178:           /* the derivative of  cnm1x * thetanm1x  wrt x variable f is 0.5 * cnm1x if f > d otherwise it is cnm1x */
1179:           /* the derivative of  cnm1  * thetanm1   wrt x variable f is 0 if f == d, otherwise it is -0.5 * cnm1 */
1180:           /* the derivative of -cnm2  * thetanm2   wrt x variable f is 0 if f == d, otherwise it is cnm2 * thetanm1 */
1181:           if (f > d) {
1182:             PetscInt f2;

1184:             p[(degidx * Nk + kidx) * npoints + pt] += mplty * 0.5 * (cnm1x - cnm1) * p[(m1idx * Nk + km1idx) * npoints + pt];
1185:             if (m2idx >= 0) {
1186:               p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm2 * thetanm1 * p[(m2idx * Nk + km1idx) * npoints + pt];
1187:               /* second derivatives of -cnm2  * thetanm2   wrt x variable f,f2 is like - 0.5 * cnm2 */
1188:               for (f2 = f; f2 < dim; f2++) {
1189:                 PetscInt km2idx, mplty2 = ktup[f2];
1190:                 PetscInt factor;

1192:                 if (!mplty2) continue;
1193:                 ktup[f2]--;
1194:                 PetscCall(PetscDTGradedOrderToIndex(dim, ktup, &km2idx));

1196:                 factor = mplty * mplty2;
1197:                 if (f == f2) factor /= 2;
1198:                 p[(degidx * Nk + kidx) * npoints + pt] -= 0.5 * factor * cnm2 * p[(m2idx * Nk + km2idx) * npoints + pt];
1199:                 ktup[f2]++;
1200:               }
1201:             }
1202:           } else {
1203:             p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm1x * p[(m1idx * Nk + km1idx) * npoints + pt];
1204:           }
1205:           ktup[f]++;
1206:         }
1207:       }
1208:     }
1209:   }
1210:   for (degidx = 0; degidx < Ndeg; degidx++) {
1211:     PetscReal scale = scales[degidx];
1212:     PetscInt  i;

1214:     for (i = 0; i < Nk * npoints; i++) p[degidx * Nk * npoints + i] *= scale;
1215:   }
1216:   PetscCall(PetscFree(scales));
1217:   PetscCall(PetscFree2(degtup, ktup));
1218:   PetscFunctionReturn(PETSC_SUCCESS);
1219: }

1221: /*@
1222:   PetscDTPTrimmedSize - The size of the trimmed polynomial space of k-forms with a given degree and form degree,
1223:   which can be evaluated in `PetscDTPTrimmedEvalJet()`.

1225:   Input Parameters:
1226: + dim - the number of variables in the multivariate polynomials
1227: . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space.
1228: - formDegree - the degree of the form

1230:   Output Parameter:
1231: - size - The number ((`dim` + `degree`) choose (`dim` + `formDegree`)) x ((`degree` + `formDegree` - 1) choose (`formDegree`))

1233:   Level: advanced

1235: .seealso: `PetscDTPTrimmedEvalJet()`
1236: @*/
1237: PetscErrorCode PetscDTPTrimmedSize(PetscInt dim, PetscInt degree, PetscInt formDegree, PetscInt *size)
1238: {
1239:   PetscInt Nrk, Nbpt; // number of trimmed polynomials

1241:   PetscFunctionBegin;
1242:   formDegree = PetscAbsInt(formDegree);
1243:   PetscCall(PetscDTBinomialInt(degree + dim, degree + formDegree, &Nbpt));
1244:   PetscCall(PetscDTBinomialInt(degree + formDegree - 1, formDegree, &Nrk));
1245:   Nbpt *= Nrk;
1246:   *size = Nbpt;
1247:   PetscFunctionReturn(PETSC_SUCCESS);
1248: }

1250: /* there was a reference implementation based on section 4.4 of Arnold, Falk & Winther (acta numerica, 2006), but it
1251:  * was inferior to this implementation */
1252: static PetscErrorCode PetscDTPTrimmedEvalJet_Internal(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[])
1253: {
1254:   PetscInt  formDegreeOrig = formDegree;
1255:   PetscBool formNegative   = (formDegreeOrig < 0) ? PETSC_TRUE : PETSC_FALSE;

1257:   PetscFunctionBegin;
1258:   formDegree = PetscAbsInt(formDegreeOrig);
1259:   if (formDegree == 0) {
1260:     PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree, jetDegree, p));
1261:     PetscFunctionReturn(PETSC_SUCCESS);
1262:   }
1263:   if (formDegree == dim) {
1264:     PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p));
1265:     PetscFunctionReturn(PETSC_SUCCESS);
1266:   }
1267:   PetscInt Nbpt;
1268:   PetscCall(PetscDTPTrimmedSize(dim, degree, formDegree, &Nbpt));
1269:   PetscInt Nf;
1270:   PetscCall(PetscDTBinomialInt(dim, formDegree, &Nf));
1271:   PetscInt Nk;
1272:   PetscCall(PetscDTBinomialInt(dim + jetDegree, dim, &Nk));
1273:   PetscCall(PetscArrayzero(p, Nbpt * Nf * Nk * npoints));

1275:   PetscInt Nbpm1; // number of scalar polynomials up to degree - 1;
1276:   PetscCall(PetscDTBinomialInt(dim + degree - 1, dim, &Nbpm1));
1277:   PetscReal *p_scalar;
1278:   PetscCall(PetscMalloc1(Nbpm1 * Nk * npoints, &p_scalar));
1279:   PetscCall(PetscDTPKDEvalJet(dim, npoints, points, degree - 1, jetDegree, p_scalar));
1280:   PetscInt total = 0;
1281:   // First add the full polynomials up to degree - 1 into the basis: take the scalar
1282:   // and copy one for each form component
1283:   for (PetscInt i = 0; i < Nbpm1; i++) {
1284:     const PetscReal *src = &p_scalar[i * Nk * npoints];
1285:     for (PetscInt f = 0; f < Nf; f++) {
1286:       PetscReal *dest = &p[(total++ * Nf + f) * Nk * npoints];
1287:       PetscCall(PetscArraycpy(dest, src, Nk * npoints));
1288:     }
1289:   }
1290:   PetscInt *form_atoms;
1291:   PetscCall(PetscMalloc1(formDegree + 1, &form_atoms));
1292:   // construct the interior product pattern
1293:   PetscInt(*pattern)[3];
1294:   PetscInt Nf1; // number of formDegree + 1 forms
1295:   PetscCall(PetscDTBinomialInt(dim, formDegree + 1, &Nf1));
1296:   PetscInt nnz = Nf1 * (formDegree + 1);
1297:   PetscCall(PetscMalloc1(Nf1 * (formDegree + 1), &pattern));
1298:   PetscCall(PetscDTAltVInteriorPattern(dim, formDegree + 1, pattern));
1299:   PetscReal centroid = (1. - dim) / (dim + 1.);
1300:   PetscInt *deriv;
1301:   PetscCall(PetscMalloc1(dim, &deriv));
1302:   for (PetscInt d = dim; d >= formDegree + 1; d--) {
1303:     PetscInt Nfd1; // number of formDegree + 1 forms in dimension d that include dx_0
1304:                    // (equal to the number of formDegree forms in dimension d-1)
1305:     PetscCall(PetscDTBinomialInt(d - 1, formDegree, &Nfd1));
1306:     // The number of homogeneous (degree-1) scalar polynomials in d variables
1307:     PetscInt Nh;
1308:     PetscCall(PetscDTBinomialInt(d - 1 + degree - 1, d - 1, &Nh));
1309:     const PetscReal *h_scalar = &p_scalar[(Nbpm1 - Nh) * Nk * npoints];
1310:     for (PetscInt b = 0; b < Nh; b++) {
1311:       const PetscReal *h_s = &h_scalar[b * Nk * npoints];
1312:       for (PetscInt f = 0; f < Nfd1; f++) {
1313:         // construct all formDegree+1 forms that start with dx_(dim - d) /\ ...
1314:         form_atoms[0] = dim - d;
1315:         PetscCall(PetscDTEnumSubset(d - 1, formDegree, f, &form_atoms[1]));
1316:         for (PetscInt i = 0; i < formDegree; i++) form_atoms[1 + i] += form_atoms[0] + 1;
1317:         PetscInt f_ind; // index of the resulting form
1318:         PetscCall(PetscDTSubsetIndex(dim, formDegree + 1, form_atoms, &f_ind));
1319:         PetscReal *p_f = &p[total++ * Nf * Nk * npoints];
1320:         for (PetscInt nz = 0; nz < nnz; nz++) {
1321:           PetscInt  i     = pattern[nz][0]; // formDegree component
1322:           PetscInt  j     = pattern[nz][1]; // (formDegree + 1) component
1323:           PetscInt  v     = pattern[nz][2]; // coordinate component
1324:           PetscReal scale = v < 0 ? -1. : 1.;

1326:           i     = formNegative ? (Nf - 1 - i) : i;
1327:           scale = (formNegative && (i & 1)) ? -scale : scale;
1328:           v     = v < 0 ? -(v + 1) : v;
1329:           if (j != f_ind) continue;
1330:           PetscReal *p_i = &p_f[i * Nk * npoints];
1331:           for (PetscInt jet = 0; jet < Nk; jet++) {
1332:             const PetscReal *h_jet = &h_s[jet * npoints];
1333:             PetscReal       *p_jet = &p_i[jet * npoints];

1335:             for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * h_jet[pt] * (points[pt * dim + v] - centroid);
1336:             PetscCall(PetscDTIndexToGradedOrder(dim, jet, deriv));
1337:             deriv[v]++;
1338:             PetscReal mult = deriv[v];
1339:             PetscInt  l;
1340:             PetscCall(PetscDTGradedOrderToIndex(dim, deriv, &l));
1341:             if (l >= Nk) continue;
1342:             p_jet = &p_i[l * npoints];
1343:             for (PetscInt pt = 0; pt < npoints; pt++) p_jet[pt] += scale * mult * h_jet[pt];
1344:             deriv[v]--;
1345:           }
1346:         }
1347:       }
1348:     }
1349:   }
1350:   PetscCheck(total == Nbpt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrectly counted P trimmed polynomials");
1351:   PetscCall(PetscFree(deriv));
1352:   PetscCall(PetscFree(pattern));
1353:   PetscCall(PetscFree(form_atoms));
1354:   PetscCall(PetscFree(p_scalar));
1355:   PetscFunctionReturn(PETSC_SUCCESS);
1356: }

1358: /*@
1359:   PetscDTPTrimmedEvalJet - Evaluate the jet (function and derivatives) of a basis of the trimmed polynomial k-forms up to
1360:   a given degree.

1362:   Input Parameters:
1363: + dim - the number of variables in the multivariate polynomials
1364: . npoints - the number of points to evaluate the polynomials at
1365: . points - [npoints x dim] array of point coordinates
1366: . degree - the degree (sum of degrees on the variables in a monomial) of the trimmed polynomial space to evaluate.
1367:            There are ((dim + degree) choose (dim + formDegree)) x ((degree + formDegree - 1) choose (formDegree)) polynomials in this space.
1368:            (You can use `PetscDTPTrimmedSize()` to compute this size.)
1369: . formDegree - the degree of the form
1370: - jetDegree - the maximum order partial derivative to evaluate in the jet.  There are ((dim + jetDegree) choose dim) partial derivatives
1371:               in the jet.  Choosing jetDegree = 0 means to evaluate just the function and no derivatives

1373:   Output Parameter:
1374: . p - an array containing the evaluations of the PKD polynomials' jets on the points.  The size is
1375:       `PetscDTPTrimmedSize()` x ((dim + formDegree) choose dim) x ((dim + k) choose dim) x npoints,
1376:       which also describes the order of the dimensions of this
1377:       four-dimensional array:
1378:         the first (slowest varying) dimension is basis function index;
1379:         the second dimension is component of the form;
1380:         the third dimension is jet index;
1381:         the fourth (fastest varying) dimension is the index of the evaluation point.

1383:   Level: advanced

1385:   Notes:
1386:   The ordering of the basis functions is not graded, so the basis functions are not nested by degree like `PetscDTPKDEvalJet()`.
1387:   The basis functions are not an L2-orthonormal basis on any particular domain.

1389:   The implementation is based on the description of the trimmed polynomials up to degree r as
1390:   the direct sum of polynomials up to degree (r-1) and the Koszul differential applied to
1391:   homogeneous polynomials of degree (r-1).

1393: .seealso: `PetscDTPKDEvalJet()`, `PetscDTPTrimmedSize()`
1394: @*/
1395: PetscErrorCode PetscDTPTrimmedEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt formDegree, PetscInt jetDegree, PetscReal p[])
1396: {
1397:   PetscFunctionBegin;
1398:   PetscCall(PetscDTPTrimmedEvalJet_Internal(dim, npoints, points, degree, formDegree, jetDegree, p));
1399:   PetscFunctionReturn(PETSC_SUCCESS);
1400: }

1402: /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V
1403:  * with lds n; diag and subdiag are overwritten */
1404: static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[], PetscReal eigs[], PetscScalar V[])
1405: {
1406:   char          jobz   = 'V'; /* eigenvalues and eigenvectors */
1407:   char          range  = 'A'; /* all eigenvalues will be found */
1408:   PetscReal     VL     = 0.;  /* ignored because range is 'A' */
1409:   PetscReal     VU     = 0.;  /* ignored because range is 'A' */
1410:   PetscBLASInt  IL     = 0;   /* ignored because range is 'A' */
1411:   PetscBLASInt  IU     = 0;   /* ignored because range is 'A' */
1412:   PetscReal     abstol = 0.;  /* unused */
1413:   PetscBLASInt  bn, bm, ldz;  /* bm will equal bn on exit */
1414:   PetscBLASInt *isuppz;
1415:   PetscBLASInt  lwork, liwork;
1416:   PetscReal     workquery;
1417:   PetscBLASInt  iworkquery;
1418:   PetscBLASInt *iwork;
1419:   PetscBLASInt  info;
1420:   PetscReal    *work = NULL;

1422:   PetscFunctionBegin;
1423: #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1424:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1425: #endif
1426:   PetscCall(PetscBLASIntCast(n, &bn));
1427:   PetscCall(PetscBLASIntCast(n, &ldz));
1428: #if !defined(PETSC_MISSING_LAPACK_STEGR)
1429:   PetscCall(PetscMalloc1(2 * n, &isuppz));
1430:   lwork  = -1;
1431:   liwork = -1;
1432:   PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, &workquery, &lwork, &iworkquery, &liwork, &info));
1433:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEGR error");
1434:   lwork  = (PetscBLASInt)workquery;
1435:   liwork = (PetscBLASInt)iworkquery;
1436:   PetscCall(PetscMalloc2(lwork, &work, liwork, &iwork));
1437:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1438:   PetscCallBLAS("LAPACKstegr", LAPACKstegr_(&jobz, &range, &bn, diag, subdiag, &VL, &VU, &IL, &IU, &abstol, &bm, eigs, V, &ldz, isuppz, work, &lwork, iwork, &liwork, &info));
1439:   PetscCall(PetscFPTrapPop());
1440:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEGR error");
1441:   PetscCall(PetscFree2(work, iwork));
1442:   PetscCall(PetscFree(isuppz));
1443: #elif !defined(PETSC_MISSING_LAPACK_STEQR)
1444:   jobz = 'I'; /* Compute eigenvalues and eigenvectors of the
1445:                  tridiagonal matrix.  Z is initialized to the identity
1446:                  matrix. */
1447:   PetscCall(PetscMalloc1(PetscMax(1, 2 * n - 2), &work));
1448:   PetscCallBLAS("LAPACKsteqr", LAPACKsteqr_("I", &bn, diag, subdiag, V, &ldz, work, &info));
1449:   PetscCall(PetscFPTrapPop());
1450:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEQR error");
1451:   PetscCall(PetscFree(work));
1452:   PetscCall(PetscArraycpy(eigs, diag, n));
1453: #endif
1454:   PetscFunctionReturn(PETSC_SUCCESS);
1455: }

1457: /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi
1458:  * quadrature rules on the interval [-1, 1] */
1459: static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw)
1460: {
1461:   PetscReal twoab1;
1462:   PetscInt  m = n - 2;
1463:   PetscReal a = alpha + 1.;
1464:   PetscReal b = beta + 1.;
1465:   PetscReal gra, grb;

1467:   PetscFunctionBegin;
1468:   twoab1 = PetscPowReal(2., a + b - 1.);
1469: #if defined(PETSC_HAVE_LGAMMA)
1470:   grb = PetscExpReal(2. * PetscLGamma(b + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + a + 1.) - (PetscLGamma(m + b + 1) + PetscLGamma(m + a + b + 1.)));
1471:   gra = PetscExpReal(2. * PetscLGamma(a + 1.) + PetscLGamma(m + 1.) + PetscLGamma(m + b + 1.) - (PetscLGamma(m + a + 1) + PetscLGamma(m + a + b + 1.)));
1472: #else
1473:   {
1474:     PetscInt alphai = (PetscInt)alpha;
1475:     PetscInt betai  = (PetscInt)beta;

1477:     if ((PetscReal)alphai == alpha && (PetscReal)betai == beta) {
1478:       PetscReal binom1, binom2;

1480:       PetscCall(PetscDTBinomial(m + b, b, &binom1));
1481:       PetscCall(PetscDTBinomial(m + a + b, b, &binom2));
1482:       grb = 1. / (binom1 * binom2);
1483:       PetscCall(PetscDTBinomial(m + a, a, &binom1));
1484:       PetscCall(PetscDTBinomial(m + a + b, a, &binom2));
1485:       gra = 1. / (binom1 * binom2);
1486:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
1487:   }
1488: #endif
1489:   *leftw  = twoab1 * grb / b;
1490:   *rightw = twoab1 * gra / a;
1491:   PetscFunctionReturn(PETSC_SUCCESS);
1492: }

1494: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
1495:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
1496: static inline PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
1497: {
1498:   PetscReal pn1, pn2;
1499:   PetscReal cnm1, cnm1x, cnm2;
1500:   PetscInt  k;

1502:   PetscFunctionBegin;
1503:   if (!n) {
1504:     *P = 1.0;
1505:     PetscFunctionReturn(PETSC_SUCCESS);
1506:   }
1507:   PetscDTJacobiRecurrence_Internal(1, a, b, cnm1, cnm1x, cnm2);
1508:   pn2 = 1.;
1509:   pn1 = cnm1 + cnm1x * x;
1510:   if (n == 1) {
1511:     *P = pn1;
1512:     PetscFunctionReturn(PETSC_SUCCESS);
1513:   }
1514:   *P = 0.0;
1515:   for (k = 2; k < n + 1; ++k) {
1516:     PetscDTJacobiRecurrence_Internal(k, a, b, cnm1, cnm1x, cnm2);

1518:     *P  = (cnm1 + cnm1x * x) * pn1 - cnm2 * pn2;
1519:     pn2 = pn1;
1520:     pn1 = *P;
1521:   }
1522:   PetscFunctionReturn(PETSC_SUCCESS);
1523: }

1525: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
1526: static inline PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P)
1527: {
1528:   PetscReal nP;
1529:   PetscInt  i;

1531:   PetscFunctionBegin;
1532:   *P = 0.0;
1533:   if (k > n) PetscFunctionReturn(PETSC_SUCCESS);
1534:   PetscCall(PetscDTComputeJacobi(a + k, b + k, n - k, x, &nP));
1535:   for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5;
1536:   *P = nP;
1537:   PetscFunctionReturn(PETSC_SUCCESS);
1538: }

1540: static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[])
1541: {
1542:   PetscInt  maxIter = 100;
1543:   PetscReal eps     = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON));
1544:   PetscReal a1, a6, gf;
1545:   PetscInt  k;

1547:   PetscFunctionBegin;

1549:   a1 = PetscPowReal(2.0, a + b + 1);
1550: #if defined(PETSC_HAVE_LGAMMA)
1551:   {
1552:     PetscReal a2, a3, a4, a5;
1553:     a2 = PetscLGamma(a + npoints + 1);
1554:     a3 = PetscLGamma(b + npoints + 1);
1555:     a4 = PetscLGamma(a + b + npoints + 1);
1556:     a5 = PetscLGamma(npoints + 1);
1557:     gf = PetscExpReal(a2 + a3 - (a4 + a5));
1558:   }
1559: #else
1560:   {
1561:     PetscInt ia, ib;

1563:     ia = (PetscInt)a;
1564:     ib = (PetscInt)b;
1565:     gf = 1.;
1566:     if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */
1567:       for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k);
1568:     } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */
1569:       for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k);
1570:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable.");
1571:   }
1572: #endif

1574:   a6 = a1 * gf;
1575:   /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
1576:    Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
1577:   for (k = 0; k < npoints; ++k) {
1578:     PetscReal r = PetscCosReal(PETSC_PI * (1. - (4. * k + 3. + 2. * b) / (4. * npoints + 2. * (a + b + 1.)))), dP;
1579:     PetscInt  j;

1581:     if (k > 0) r = 0.5 * (r + x[k - 1]);
1582:     for (j = 0; j < maxIter; ++j) {
1583:       PetscReal s = 0.0, delta, f, fp;
1584:       PetscInt  i;

1586:       for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
1587:       PetscCall(PetscDTComputeJacobi(a, b, npoints, r, &f));
1588:       PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp));
1589:       delta = f / (fp - f * s);
1590:       r     = r - delta;
1591:       if (PetscAbsReal(delta) < eps) break;
1592:     }
1593:     x[k] = r;
1594:     PetscCall(PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP));
1595:     w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
1596:   }
1597:   PetscFunctionReturn(PETSC_SUCCESS);
1598: }

1600: /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi
1601:  * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */
1602: static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s)
1603: {
1604:   PetscInt i;

1606:   PetscFunctionBegin;
1607:   for (i = 0; i < nPoints; i++) {
1608:     PetscReal A, B, C;

1610:     PetscDTJacobiRecurrence_Internal(i + 1, a, b, A, B, C);
1611:     d[i] = -A / B;
1612:     if (i) s[i - 1] *= C / B;
1613:     if (i < nPoints - 1) s[i] = 1. / B;
1614:   }
1615:   PetscFunctionReturn(PETSC_SUCCESS);
1616: }

1618: static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
1619: {
1620:   PetscReal mu0;
1621:   PetscReal ga, gb, gab;
1622:   PetscInt  i;

1624:   PetscFunctionBegin;
1625:   PetscCall(PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite));

1627: #if defined(PETSC_HAVE_TGAMMA)
1628:   ga  = PetscTGamma(a + 1);
1629:   gb  = PetscTGamma(b + 1);
1630:   gab = PetscTGamma(a + b + 2);
1631: #else
1632:   {
1633:     PetscInt ia, ib;

1635:     ia = (PetscInt)a;
1636:     ib = (PetscInt)b;
1637:     if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */
1638:       PetscCall(PetscDTFactorial(ia, &ga));
1639:       PetscCall(PetscDTFactorial(ib, &gb));
1640:       PetscCall(PetscDTFactorial(ia + ib + 1, &gab));
1641:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "tgamma() - math routine is unavailable.");
1642:   }
1643: #endif
1644:   mu0 = PetscPowReal(2., a + b + 1.) * ga * gb / gab;

1646: #if defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1647:   {
1648:     PetscReal   *diag, *subdiag;
1649:     PetscScalar *V;

1651:     PetscCall(PetscMalloc2(npoints, &diag, npoints, &subdiag));
1652:     PetscCall(PetscMalloc1(npoints * npoints, &V));
1653:     PetscCall(PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag));
1654:     for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]);
1655:     PetscCall(PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V));
1656:     for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0;
1657:     PetscCall(PetscFree(V));
1658:     PetscCall(PetscFree2(diag, subdiag));
1659:   }
1660: #else
1661:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1662: #endif
1663:   { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the
1664:        eigenvalues are not guaranteed to be in ascending order.  So we heave a passive aggressive sigh and check that
1665:        the eigenvalues are sorted */
1666:     PetscBool sorted;

1668:     PetscCall(PetscSortedReal(npoints, x, &sorted));
1669:     if (!sorted) {
1670:       PetscInt  *order, i;
1671:       PetscReal *tmp;

1673:       PetscCall(PetscMalloc2(npoints, &order, npoints, &tmp));
1674:       for (i = 0; i < npoints; i++) order[i] = i;
1675:       PetscCall(PetscSortRealWithPermutation(npoints, x, order));
1676:       PetscCall(PetscArraycpy(tmp, x, npoints));
1677:       for (i = 0; i < npoints; i++) x[i] = tmp[order[i]];
1678:       PetscCall(PetscArraycpy(tmp, w, npoints));
1679:       for (i = 0; i < npoints; i++) w[i] = tmp[order[i]];
1680:       PetscCall(PetscFree2(order, tmp));
1681:     }
1682:   }
1683:   PetscFunctionReturn(PETSC_SUCCESS);
1684: }

1686: static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1687: {
1688:   PetscFunctionBegin;
1689:   PetscCheck(npoints >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive");
1690:   /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1691:   PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1.");
1692:   PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1.");

1694:   if (newton) PetscCall(PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w));
1695:   else PetscCall(PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w));
1696:   if (alpha == beta) { /* symmetrize */
1697:     PetscInt i;
1698:     for (i = 0; i < (npoints + 1) / 2; i++) {
1699:       PetscInt  j  = npoints - 1 - i;
1700:       PetscReal xi = x[i];
1701:       PetscReal xj = x[j];
1702:       PetscReal wi = w[i];
1703:       PetscReal wj = w[j];

1705:       x[i] = (xi - xj) / 2.;
1706:       x[j] = (xj - xi) / 2.;
1707:       w[i] = w[j] = (wi + wj) / 2.;
1708:     }
1709:   }
1710:   PetscFunctionReturn(PETSC_SUCCESS);
1711: }

1713: /*@
1714:   PetscDTGaussJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1715:   $(x-a)^\alpha (x-b)^\beta$.

1717:   Not Collective

1719:   Input Parameters:
1720: + npoints - the number of points in the quadrature rule
1721: . a - the left endpoint of the interval
1722: . b - the right endpoint of the interval
1723: . alpha - the left exponent
1724: - beta - the right exponent

1726:   Output Parameters:
1727: + x - array of length `npoints`, the locations of the quadrature points
1728: - w - array of length `npoints`, the weights of the quadrature points

1730:   Level: intermediate

1732:   Note:
1733:   This quadrature rule is exact for polynomials up to degree 2*npoints - 1.

1735: .seealso: `PetscDTGaussQuadrature()`
1736: @*/
1737: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1738: {
1739:   PetscInt i;

1741:   PetscFunctionBegin;
1742:   PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal));
1743:   if (a != -1. || b != 1.) { /* shift */
1744:     for (i = 0; i < npoints; i++) {
1745:       x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1746:       w[i] *= (b - a) / 2.;
1747:     }
1748:   }
1749:   PetscFunctionReturn(PETSC_SUCCESS);
1750: }

1752: static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1753: {
1754:   PetscInt i;

1756:   PetscFunctionBegin;
1757:   PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be positive");
1758:   /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1759:   PetscCheck(alpha > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "alpha must be > -1.");
1760:   PetscCheck(beta > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "beta must be > -1.");

1762:   x[0]           = -1.;
1763:   x[npoints - 1] = 1.;
1764:   if (npoints > 2) PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints - 2, alpha + 1., beta + 1., &x[1], &w[1], newton));
1765:   for (i = 1; i < npoints - 1; i++) w[i] /= (1. - x[i] * x[i]);
1766:   PetscCall(PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints - 1]));
1767:   PetscFunctionReturn(PETSC_SUCCESS);
1768: }

1770: /*@
1771:   PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1772:   $(x-a)^\alpha (x-b)^\beta$, with endpoints a and b included as quadrature points.

1774:   Not Collective

1776:   Input Parameters:
1777: + npoints - the number of points in the quadrature rule
1778: . a - the left endpoint of the interval
1779: . b - the right endpoint of the interval
1780: . alpha - the left exponent
1781: - beta - the right exponent

1783:   Output Parameters:
1784: + x - array of length `npoints`, the locations of the quadrature points
1785: - w - array of length `npoints`, the weights of the quadrature points

1787:   Level: intermediate

1789:   Note:
1790:   This quadrature rule is exact for polynomials up to degree 2*npoints - 3.

1792: .seealso: `PetscDTGaussJacobiQuadrature()`
1793: @*/
1794: PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1795: {
1796:   PetscInt i;

1798:   PetscFunctionBegin;
1799:   PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal));
1800:   if (a != -1. || b != 1.) { /* shift */
1801:     for (i = 0; i < npoints; i++) {
1802:       x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1803:       w[i] *= (b - a) / 2.;
1804:     }
1805:   }
1806:   PetscFunctionReturn(PETSC_SUCCESS);
1807: }

1809: /*@
1810:    PetscDTGaussQuadrature - create Gauss-Legendre quadrature

1812:    Not Collective

1814:    Input Parameters:
1815: +  npoints - number of points
1816: .  a - left end of interval (often-1)
1817: -  b - right end of interval (often +1)

1819:    Output Parameters:
1820: +  x - quadrature points
1821: -  w - quadrature weights

1823:    Level: intermediate

1825:    References:
1826: .  * - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.

1828: .seealso: `PetscDTLegendreEval()`, `PetscDTGaussJacobiQuadrature()`
1829: @*/
1830: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
1831: {
1832:   PetscInt i;

1834:   PetscFunctionBegin;
1835:   PetscCall(PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal));
1836:   if (a != -1. || b != 1.) { /* shift */
1837:     for (i = 0; i < npoints; i++) {
1838:       x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1839:       w[i] *= (b - a) / 2.;
1840:     }
1841:   }
1842:   PetscFunctionReturn(PETSC_SUCCESS);
1843: }

1845: /*@C
1846:    PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre
1847:                       nodes of a given size on the domain [-1,1]

1849:    Not Collective

1851:    Input Parameters:
1852: +  n - number of grid nodes
1853: -  type - `PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA` or `PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON`

1855:    Output Parameters:
1856: +  x - quadrature points
1857: -  w - quadrature weights

1859:    Level: intermediate

1861:    Notes:
1862:     For n > 30  the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
1863:           close enough to the desired solution

1865:    These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes

1867:    See  https://epubs.siam.org/doi/abs/10.1137/110855442  https://epubs.siam.org/doi/abs/10.1137/120889873 for better ways to compute GLL nodes

1869: .seealso: `PetscDTGaussQuadrature()`, `PetscGaussLobattoLegendreCreateType`

1871: @*/
1872: PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints, PetscGaussLobattoLegendreCreateType type, PetscReal *x, PetscReal *w)
1873: {
1874:   PetscBool newton;

1876:   PetscFunctionBegin;
1877:   PetscCheck(npoints >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide at least 2 grid points per element");
1878:   newton = (PetscBool)(type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON);
1879:   PetscCall(PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton));
1880:   PetscFunctionReturn(PETSC_SUCCESS);
1881: }

1883: /*@
1884:   PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature

1886:   Not Collective

1888:   Input Parameters:
1889: + dim     - The spatial dimension
1890: . Nc      - The number of components
1891: . npoints - number of points in one dimension
1892: . a       - left end of interval (often-1)
1893: - b       - right end of interval (often +1)

1895:   Output Parameter:
1896: . q - A `PetscQuadrature` object

1898:   Level: intermediate

1900: .seealso: `PetscDTGaussQuadrature()`, `PetscDTLegendreEval()`
1901: @*/
1902: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1903: {
1904:   DMPolytopeType ct;
1905:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints * PetscSqr(npoints) : PetscSqr(npoints) : npoints;
1906:   PetscReal     *x, *w, *xw, *ww;

1908:   PetscFunctionBegin;
1909:   PetscCall(PetscMalloc1(totpoints * dim, &x));
1910:   PetscCall(PetscMalloc1(totpoints * Nc, &w));
1911:   /* Set up the Golub-Welsch system */
1912:   switch (dim) {
1913:   case 0:
1914:     ct = DM_POLYTOPE_POINT;
1915:     PetscCall(PetscFree(x));
1916:     PetscCall(PetscFree(w));
1917:     PetscCall(PetscMalloc1(1, &x));
1918:     PetscCall(PetscMalloc1(Nc, &w));
1919:     x[0] = 0.0;
1920:     for (PetscInt c = 0; c < Nc; ++c) w[c] = 1.0;
1921:     break;
1922:   case 1:
1923:     ct = DM_POLYTOPE_SEGMENT;
1924:     PetscCall(PetscMalloc1(npoints, &ww));
1925:     PetscCall(PetscDTGaussQuadrature(npoints, a, b, x, ww));
1926:     for (PetscInt i = 0; i < npoints; ++i)
1927:       for (PetscInt c = 0; c < Nc; ++c) w[i * Nc + c] = ww[i];
1928:     PetscCall(PetscFree(ww));
1929:     break;
1930:   case 2:
1931:     ct = DM_POLYTOPE_QUADRILATERAL;
1932:     PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww));
1933:     PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww));
1934:     for (PetscInt i = 0; i < npoints; ++i) {
1935:       for (PetscInt j = 0; j < npoints; ++j) {
1936:         x[(i * npoints + j) * dim + 0] = xw[i];
1937:         x[(i * npoints + j) * dim + 1] = xw[j];
1938:         for (PetscInt c = 0; c < Nc; ++c) w[(i * npoints + j) * Nc + c] = ww[i] * ww[j];
1939:       }
1940:     }
1941:     PetscCall(PetscFree2(xw, ww));
1942:     break;
1943:   case 3:
1944:     ct = DM_POLYTOPE_HEXAHEDRON;
1945:     PetscCall(PetscMalloc2(npoints, &xw, npoints, &ww));
1946:     PetscCall(PetscDTGaussQuadrature(npoints, a, b, xw, ww));
1947:     for (PetscInt i = 0; i < npoints; ++i) {
1948:       for (PetscInt j = 0; j < npoints; ++j) {
1949:         for (PetscInt k = 0; k < npoints; ++k) {
1950:           x[((i * npoints + j) * npoints + k) * dim + 0] = xw[i];
1951:           x[((i * npoints + j) * npoints + k) * dim + 1] = xw[j];
1952:           x[((i * npoints + j) * npoints + k) * dim + 2] = xw[k];
1953:           for (PetscInt c = 0; c < Nc; ++c) w[((i * npoints + j) * npoints + k) * Nc + c] = ww[i] * ww[j] * ww[k];
1954:         }
1955:       }
1956:     }
1957:     PetscCall(PetscFree2(xw, ww));
1958:     break;
1959:   default:
1960:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %" PetscInt_FMT, dim);
1961:   }
1962:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
1963:   PetscCall(PetscQuadratureSetCellType(*q, ct));
1964:   PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1));
1965:   PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w));
1966:   PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "GaussTensor"));
1967:   PetscFunctionReturn(PETSC_SUCCESS);
1968: }

1970: /*@
1971:   PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex

1973:   Not Collective

1975:   Input Parameters:
1976: + dim     - The simplex dimension
1977: . Nc      - The number of components
1978: . npoints - The number of points in one dimension
1979: . a       - left end of interval (often-1)
1980: - b       - right end of interval (often +1)

1982:   Output Parameter:
1983: . q - A `PetscQuadrature` object

1985:   Level: intermediate

1987:   Note:
1988:   For `dim` == 1, this is Gauss-Legendre quadrature

1990:   References:
1991: . * - Karniadakis and Sherwin.  FIAT

1993: .seealso: `PetscDTGaussTensorQuadrature()`, `PetscDTGaussQuadrature()`
1994: @*/
1995: PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1996: {
1997:   DMPolytopeType ct;
1998:   PetscInt       totpoints;
1999:   PetscReal     *p1, *w1;
2000:   PetscReal     *x, *w;

2002:   PetscFunctionBegin;
2003:   PetscCheck(!(a != -1.0) && !(b != 1.0), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
2004:   switch (dim) {
2005:   case 0:
2006:     ct = DM_POLYTOPE_POINT;
2007:     break;
2008:   case 1:
2009:     ct = DM_POLYTOPE_SEGMENT;
2010:     break;
2011:   case 2:
2012:     ct = DM_POLYTOPE_TRIANGLE;
2013:     break;
2014:   case 3:
2015:     ct = DM_POLYTOPE_TETRAHEDRON;
2016:     break;
2017:   default:
2018:     ct = DM_POLYTOPE_UNKNOWN;
2019:   }
2020:   totpoints = 1;
2021:   for (PetscInt i = 0; i < dim; ++i) totpoints *= npoints;
2022:   PetscCall(PetscMalloc1(totpoints * dim, &x));
2023:   PetscCall(PetscMalloc1(totpoints * Nc, &w));
2024:   PetscCall(PetscMalloc2(npoints, &p1, npoints, &w1));
2025:   for (PetscInt i = 0; i < totpoints * Nc; ++i) w[i] = 1.;
2026:   for (PetscInt i = 0, totprev = 1, totrem = totpoints / npoints; i < dim; ++i) {
2027:     PetscReal mul;

2029:     mul = PetscPowReal(2., -i);
2030:     PetscCall(PetscDTGaussJacobiQuadrature(npoints, -1., 1., i, 0.0, p1, w1));
2031:     for (PetscInt pt = 0, l = 0; l < totprev; l++) {
2032:       for (PetscInt j = 0; j < npoints; j++) {
2033:         for (PetscInt m = 0; m < totrem; m++, pt++) {
2034:           for (PetscInt k = 0; k < i; k++) x[pt * dim + k] = (x[pt * dim + k] + 1.) * (1. - p1[j]) * 0.5 - 1.;
2035:           x[pt * dim + i] = p1[j];
2036:           for (PetscInt c = 0; c < Nc; c++) w[pt * Nc + c] *= mul * w1[j];
2037:         }
2038:       }
2039:     }
2040:     totprev *= npoints;
2041:     totrem /= npoints;
2042:   }
2043:   PetscCall(PetscFree2(p1, w1));
2044:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
2045:   PetscCall(PetscQuadratureSetCellType(*q, ct));
2046:   PetscCall(PetscQuadratureSetOrder(*q, 2 * npoints - 1));
2047:   PetscCall(PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w));
2048:   PetscCall(PetscObjectChangeTypeName((PetscObject)*q, "StroudConical"));
2049:   PetscFunctionReturn(PETSC_SUCCESS);
2050: }

2052: static PetscBool MinSymTriQuadCite       = PETSC_FALSE;
2053: const char       MinSymTriQuadCitation[] = "@article{WitherdenVincent2015,\n"
2054:                                            "  title = {On the identification of symmetric quadrature rules for finite element methods},\n"
2055:                                            "  journal = {Computers & Mathematics with Applications},\n"
2056:                                            "  volume = {69},\n"
2057:                                            "  number = {10},\n"
2058:                                            "  pages = {1232-1241},\n"
2059:                                            "  year = {2015},\n"
2060:                                            "  issn = {0898-1221},\n"
2061:                                            "  doi = {10.1016/j.camwa.2015.03.017},\n"
2062:                                            "  url = {https://www.sciencedirect.com/science/article/pii/S0898122115001224},\n"
2063:                                            "  author = {F.D. Witherden and P.E. Vincent},\n"
2064:                                            "}\n";

2066: #include "petscdttriquadrules.h"

2068: static PetscBool MinSymTetQuadCite       = PETSC_FALSE;
2069: const char       MinSymTetQuadCitation[] = "@article{JaskowiecSukumar2021\n"
2070:                                            "  author = {Jaskowiec, Jan and Sukumar, N.},\n"
2071:                                            "  title = {High-order symmetric cubature rules for tetrahedra and pyramids},\n"
2072:                                            "  journal = {International Journal for Numerical Methods in Engineering},\n"
2073:                                            "  volume = {122},\n"
2074:                                            "  number = {1},\n"
2075:                                            "  pages = {148-171},\n"
2076:                                            "  doi = {10.1002/nme.6528},\n"
2077:                                            "  url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6528},\n"
2078:                                            "  eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6528},\n"
2079:                                            "  year = {2021}\n"
2080:                                            "}\n";

2082: #include "petscdttetquadrules.h"

2084: // https://en.wikipedia.org/wiki/Partition_(number_theory)
2085: static PetscErrorCode PetscDTPartitionNumber(PetscInt n, PetscInt *p)
2086: {
2087:   // sequence A000041 in the OEIS
2088:   const PetscInt partition[]   = {1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135, 176, 231, 297, 385, 490, 627, 792, 1002, 1255, 1575, 1958, 2436, 3010, 3718, 4565, 5604};
2089:   PetscInt       tabulated_max = PETSC_STATIC_ARRAY_LENGTH(partition) - 1;

2091:   PetscFunctionBegin;
2092:   PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Partition number not defined for negative number %" PetscInt_FMT, n);
2093:   // not implementing the pentagonal number recurrence, we don't need partition numbers for n that high
2094:   PetscCheck(n <= tabulated_max, PETSC_COMM_SELF, PETSC_ERR_SUP, "Partition numbers only tabulated up to %" PetscInt_FMT ", not computed for %" PetscInt_FMT, tabulated_max, n);
2095:   *p = partition[n];
2096:   PetscFunctionReturn(PETSC_SUCCESS);
2097: }

2099: /*@
2100:   PetscDTSimplexQuadrature - Create a quadrature rule for a simplex that exactly integrates polynomials up to a given degree.

2102:   Not Collective

2104:   Input Parameters:
2105: + dim     - The spatial dimension of the simplex (1 = segment, 2 = triangle, 3 = tetrahedron)
2106: . degree  - The largest polynomial degree that is required to be integrated exactly
2107: - type    - left end of interval (often-1)

2109:   Output Parameter:
2110: . quad    - A `PetscQuadrature` object for integration over the biunit simplex
2111:             (defined by the bounds $x_i >= -1$ and $\sum_i x_i <= 2 - d$) that is exact for
2112:             polynomials up to the given degree

2114:   Level: intermediate

2116: .seealso: `PetscDTSimplexQuadratureType`, `PetscDTGaussQuadrature()`, `PetscDTStroudCononicalQuadrature()`, `PetscQuadrature`
2117: @*/
2118: PetscErrorCode PetscDTSimplexQuadrature(PetscInt dim, PetscInt degree, PetscDTSimplexQuadratureType type, PetscQuadrature *quad)
2119: {
2120:   PetscDTSimplexQuadratureType orig_type = type;

2122:   PetscFunctionBegin;
2123:   PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative dimension %" PetscInt_FMT, dim);
2124:   PetscCheck(degree >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %" PetscInt_FMT, degree);
2125:   if (type == PETSCDTSIMPLEXQUAD_DEFAULT) type = PETSCDTSIMPLEXQUAD_MINSYM;
2126:   if (type == PETSCDTSIMPLEXQUAD_CONIC || dim < 2) {
2127:     PetscInt points_per_dim = (degree + 2) / 2; // ceil((degree + 1) / 2);
2128:     PetscCall(PetscDTStroudConicalQuadrature(dim, 1, points_per_dim, -1, 1, quad));
2129:   } else {
2130:     DMPolytopeType    ct;
2131:     PetscInt          n    = dim + 1;
2132:     PetscInt          fact = 1;
2133:     PetscInt         *part, *perm;
2134:     PetscInt          p = 0;
2135:     PetscInt          max_degree;
2136:     const PetscInt   *nodes_per_type     = NULL;
2137:     const PetscInt   *all_num_full_nodes = NULL;
2138:     const PetscReal **weights_list       = NULL;
2139:     const PetscReal **compact_nodes_list = NULL;
2140:     const char       *citation           = NULL;
2141:     PetscBool        *cited              = NULL;

2143:     switch (dim) {
2144:     case 0:
2145:       ct = DM_POLYTOPE_POINT;
2146:       break;
2147:     case 1:
2148:       ct = DM_POLYTOPE_SEGMENT;
2149:       break;
2150:     case 2:
2151:       ct = DM_POLYTOPE_TRIANGLE;
2152:       break;
2153:     case 3:
2154:       ct = DM_POLYTOPE_TETRAHEDRON;
2155:       break;
2156:     default:
2157:       ct = DM_POLYTOPE_UNKNOWN;
2158:     }
2159:     switch (dim) {
2160:     case 2:
2161:       cited              = &MinSymTriQuadCite;
2162:       citation           = MinSymTriQuadCitation;
2163:       max_degree         = PetscDTWVTriQuad_max_degree;
2164:       nodes_per_type     = PetscDTWVTriQuad_num_orbits;
2165:       all_num_full_nodes = PetscDTWVTriQuad_num_nodes;
2166:       weights_list       = PetscDTWVTriQuad_weights;
2167:       compact_nodes_list = PetscDTWVTriQuad_orbits;
2168:       break;
2169:     case 3:
2170:       cited              = &MinSymTetQuadCite;
2171:       citation           = MinSymTetQuadCitation;
2172:       max_degree         = PetscDTJSTetQuad_max_degree;
2173:       nodes_per_type     = PetscDTJSTetQuad_num_orbits;
2174:       all_num_full_nodes = PetscDTJSTetQuad_num_nodes;
2175:       weights_list       = PetscDTJSTetQuad_weights;
2176:       compact_nodes_list = PetscDTJSTetQuad_orbits;
2177:       break;
2178:     default:
2179:       max_degree = -1;
2180:       break;
2181:     }

2183:     if (degree > max_degree) {
2184:       if (orig_type == PETSCDTSIMPLEXQUAD_DEFAULT) {
2185:         // fall back to conic
2186:         PetscCall(PetscDTSimplexQuadrature(dim, degree, PETSCDTSIMPLEXQUAD_CONIC, quad));
2187:         PetscFunctionReturn(PETSC_SUCCESS);
2188:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Minimal symmetric quadrature for dim %" PetscInt_FMT ", degree %" PetscInt_FMT " unsupported", dim, degree);
2189:     }

2191:     PetscCall(PetscCitationsRegister(citation, cited));

2193:     PetscCall(PetscDTPartitionNumber(n, &p));
2194:     for (PetscInt d = 2; d <= n; d++) fact *= d;

2196:     PetscInt         num_full_nodes      = all_num_full_nodes[degree];
2197:     const PetscReal *all_compact_nodes   = compact_nodes_list[degree];
2198:     const PetscReal *all_compact_weights = weights_list[degree];
2199:     nodes_per_type                       = &nodes_per_type[p * degree];

2201:     PetscReal      *points;
2202:     PetscReal      *counts;
2203:     PetscReal      *weights;
2204:     PetscReal      *bary_to_biunit; // row-major transformation of barycentric coordinate to biunit
2205:     PetscQuadrature q;

2207:     // compute the transformation
2208:     PetscCall(PetscMalloc1(n * dim, &bary_to_biunit));
2209:     for (PetscInt d = 0; d < dim; d++) {
2210:       for (PetscInt b = 0; b < n; b++) bary_to_biunit[d * n + b] = (d == b) ? 1.0 : -1.0;
2211:     }

2213:     PetscCall(PetscMalloc3(n, &part, n, &perm, n, &counts));
2214:     PetscCall(PetscCalloc1(num_full_nodes * dim, &points));
2215:     PetscCall(PetscMalloc1(num_full_nodes, &weights));

2217:     // (0, 0, ...) is the first partition lexicographically
2218:     PetscCall(PetscArrayzero(part, n));
2219:     PetscCall(PetscArrayzero(counts, n));
2220:     counts[0] = n;

2222:     // for each partition
2223:     for (PetscInt s = 0, node_offset = 0; s < p; s++) {
2224:       PetscInt num_compact_coords = part[n - 1] + 1;

2226:       const PetscReal *compact_nodes   = all_compact_nodes;
2227:       const PetscReal *compact_weights = all_compact_weights;
2228:       all_compact_nodes += num_compact_coords * nodes_per_type[s];
2229:       all_compact_weights += nodes_per_type[s];

2231:       // for every permutation of the vertices
2232:       for (PetscInt f = 0; f < fact; f++) {
2233:         PetscCall(PetscDTEnumPerm(n, f, perm, NULL));

2235:         // check if it is a valid permutation
2236:         PetscInt digit;
2237:         for (digit = 1; digit < n; digit++) {
2238:           // skip permutations that would duplicate a node because it has a smaller symmetry group
2239:           if (part[digit - 1] == part[digit] && perm[digit - 1] > perm[digit]) break;
2240:         }
2241:         if (digit < n) continue;

2243:         // create full nodes from this permutation of the compact nodes
2244:         PetscReal *full_nodes   = &points[node_offset * dim];
2245:         PetscReal *full_weights = &weights[node_offset];

2247:         PetscCall(PetscArraycpy(full_weights, compact_weights, nodes_per_type[s]));
2248:         for (PetscInt b = 0; b < n; b++) {
2249:           for (PetscInt d = 0; d < dim; d++) {
2250:             for (PetscInt node = 0; node < nodes_per_type[s]; node++) full_nodes[node * dim + d] += bary_to_biunit[d * n + perm[b]] * compact_nodes[node * num_compact_coords + part[b]];
2251:           }
2252:         }
2253:         node_offset += nodes_per_type[s];
2254:       }

2256:       if (s < p - 1) { // Generate the next partition
2257:         /* A partition is described by the number of coordinates that are in
2258:          * each set of duplicates (counts) and redundantly by mapping each
2259:          * index to its set of duplicates (part)
2260:          *
2261:          * Counts should always be in nonincreasing order
2262:          *
2263:          * We want to generate the partitions lexically by part, which means
2264:          * finding the last index where count > 1 and reducing by 1.
2265:          *
2266:          * For the new counts beyond that index, we eagerly assign the remaining
2267:          * capacity of the partition to smaller indices (ensures lexical ordering),
2268:          * while respecting the nonincreasing invariant of the counts
2269:          */
2270:         PetscInt last_digit            = part[n - 1];
2271:         PetscInt last_digit_with_extra = last_digit;
2272:         while (counts[last_digit_with_extra] == 1) last_digit_with_extra--;
2273:         PetscInt limit               = --counts[last_digit_with_extra];
2274:         PetscInt total_to_distribute = last_digit - last_digit_with_extra + 1;
2275:         for (PetscInt digit = last_digit_with_extra + 1; digit < n; digit++) {
2276:           counts[digit] = PetscMin(limit, total_to_distribute);
2277:           total_to_distribute -= PetscMin(limit, total_to_distribute);
2278:         }
2279:         for (PetscInt digit = 0, offset = 0; digit < n; digit++) {
2280:           PetscInt count = counts[digit];
2281:           for (PetscInt c = 0; c < count; c++) part[offset++] = digit;
2282:         }
2283:       }
2284:     }
2285:     PetscCall(PetscFree3(part, perm, counts));
2286:     PetscCall(PetscFree(bary_to_biunit));
2287:     PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &q));
2288:     PetscCall(PetscQuadratureSetCellType(q, ct));
2289:     PetscCall(PetscQuadratureSetOrder(q, degree));
2290:     PetscCall(PetscQuadratureSetData(q, dim, 1, num_full_nodes, points, weights));
2291:     *quad = q;
2292:   }
2293:   PetscFunctionReturn(PETSC_SUCCESS);
2294: }

2296: /*@
2297:   PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell

2299:   Not Collective

2301:   Input Parameters:
2302: + dim   - The cell dimension
2303: . level - The number of points in one dimension, 2^l
2304: . a     - left end of interval (often-1)
2305: - b     - right end of interval (often +1)

2307:   Output Parameter:
2308: . q - A `PetscQuadrature` object

2310:   Level: intermediate

2312: .seealso: `PetscDTGaussTensorQuadrature()`, `PetscQuadrature`
2313: @*/
2314: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
2315: {
2316:   DMPolytopeType  ct;
2317:   const PetscInt  p     = 16;                        /* Digits of precision in the evaluation */
2318:   const PetscReal alpha = (b - a) / 2.;              /* Half-width of the integration interval */
2319:   const PetscReal beta  = (b + a) / 2.;              /* Center of the integration interval */
2320:   const PetscReal h     = PetscPowReal(2.0, -level); /* Step size, length between x_k */
2321:   PetscReal       xk;                                /* Quadrature point x_k on reference domain [-1, 1] */
2322:   PetscReal       wk = 0.5 * PETSC_PI;               /* Quadrature weight at x_k */
2323:   PetscReal      *x, *w;
2324:   PetscInt        K, k, npoints;

2326:   PetscFunctionBegin;
2327:   PetscCheck(dim <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not yet implemented", dim);
2328:   PetscCheck(level, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
2329:   switch (dim) {
2330:   case 0:
2331:     ct = DM_POLYTOPE_POINT;
2332:     break;
2333:   case 1:
2334:     ct = DM_POLYTOPE_SEGMENT;
2335:     break;
2336:   case 2:
2337:     ct = DM_POLYTOPE_QUADRILATERAL;
2338:     break;
2339:   case 3:
2340:     ct = DM_POLYTOPE_HEXAHEDRON;
2341:     break;
2342:   default:
2343:     ct = DM_POLYTOPE_UNKNOWN;
2344:   }
2345:   /* Find K such that the weights are < 32 digits of precision */
2346:   for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2 * p; ++K) wk = 0.5 * h * PETSC_PI * PetscCoshReal(K * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(K * h)));
2347:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
2348:   PetscCall(PetscQuadratureSetCellType(*q, ct));
2349:   PetscCall(PetscQuadratureSetOrder(*q, 2 * K + 1));
2350:   npoints = 2 * K - 1;
2351:   PetscCall(PetscMalloc1(npoints * dim, &x));
2352:   PetscCall(PetscMalloc1(npoints, &w));
2353:   /* Center term */
2354:   x[0] = beta;
2355:   w[0] = 0.5 * alpha * PETSC_PI;
2356:   for (k = 1; k < K; ++k) {
2357:     wk           = 0.5 * alpha * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2358:     xk           = PetscTanhReal(0.5 * PETSC_PI * PetscSinhReal(k * h));
2359:     x[2 * k - 1] = -alpha * xk + beta;
2360:     w[2 * k - 1] = wk;
2361:     x[2 * k + 0] = alpha * xk + beta;
2362:     w[2 * k + 0] = wk;
2363:   }
2364:   PetscCall(PetscQuadratureSetData(*q, dim, 1, npoints, x, w));
2365:   PetscFunctionReturn(PETSC_SUCCESS);
2366: }

2368: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2369: {
2370:   const PetscInt  p     = 16;           /* Digits of precision in the evaluation */
2371:   const PetscReal alpha = (b - a) / 2.; /* Half-width of the integration interval */
2372:   const PetscReal beta  = (b + a) / 2.; /* Center of the integration interval */
2373:   PetscReal       h     = 1.0;          /* Step size, length between x_k */
2374:   PetscInt        l     = 0;            /* Level of refinement, h = 2^{-l} */
2375:   PetscReal       osum  = 0.0;          /* Integral on last level */
2376:   PetscReal       psum  = 0.0;          /* Integral on the level before the last level */
2377:   PetscReal       sum;                  /* Integral on current level */
2378:   PetscReal       yk;                   /* Quadrature point 1 - x_k on reference domain [-1, 1] */
2379:   PetscReal       lx, rx;               /* Quadrature points to the left and right of 0 on the real domain [a, b] */
2380:   PetscReal       wk;                   /* Quadrature weight at x_k */
2381:   PetscReal       lval, rval;           /* Terms in the quadature sum to the left and right of 0 */
2382:   PetscInt        d;                    /* Digits of precision in the integral */

2384:   PetscFunctionBegin;
2385:   PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
2386:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2387:   /* Center term */
2388:   func(&beta, ctx, &lval);
2389:   sum = 0.5 * alpha * PETSC_PI * lval;
2390:   /* */
2391:   do {
2392:     PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
2393:     PetscInt  k = 1;

2395:     ++l;
2396:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */
2397:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
2398:     psum = osum;
2399:     osum = sum;
2400:     h *= 0.5;
2401:     sum *= 0.5;
2402:     do {
2403:       wk = 0.5 * h * PETSC_PI * PetscCoshReal(k * h) / PetscSqr(PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2404:       yk = 1.0 / (PetscExpReal(0.5 * PETSC_PI * PetscSinhReal(k * h)) * PetscCoshReal(0.5 * PETSC_PI * PetscSinhReal(k * h)));
2405:       lx = -alpha * (1.0 - yk) + beta;
2406:       rx = alpha * (1.0 - yk) + beta;
2407:       func(&lx, ctx, &lval);
2408:       func(&rx, ctx, &rval);
2409:       lterm   = alpha * wk * lval;
2410:       maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
2411:       sum += lterm;
2412:       rterm   = alpha * wk * rval;
2413:       maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
2414:       sum += rterm;
2415:       ++k;
2416:       /* Only need to evaluate every other point on refined levels */
2417:       if (l != 1) ++k;
2418:     } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */

2420:     d1 = PetscLog10Real(PetscAbsReal(sum - osum));
2421:     d2 = PetscLog10Real(PetscAbsReal(sum - psum));
2422:     d3 = PetscLog10Real(maxTerm) - p;
2423:     if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
2424:     else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
2425:     d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4)));
2426:   } while (d < digits && l < 12);
2427:   *sol = sum;
2428:   PetscCall(PetscFPTrapPop());
2429:   PetscFunctionReturn(PETSC_SUCCESS);
2430: }

2432: #if defined(PETSC_HAVE_MPFR)
2433: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2434: {
2435:   const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */
2436:   PetscInt       l            = 0; /* Level of refinement, h = 2^{-l} */
2437:   mpfr_t         alpha;            /* Half-width of the integration interval */
2438:   mpfr_t         beta;             /* Center of the integration interval */
2439:   mpfr_t         h;                /* Step size, length between x_k */
2440:   mpfr_t         osum;             /* Integral on last level */
2441:   mpfr_t         psum;             /* Integral on the level before the last level */
2442:   mpfr_t         sum;              /* Integral on current level */
2443:   mpfr_t         yk;               /* Quadrature point 1 - x_k on reference domain [-1, 1] */
2444:   mpfr_t         lx, rx;           /* Quadrature points to the left and right of 0 on the real domain [a, b] */
2445:   mpfr_t         wk;               /* Quadrature weight at x_k */
2446:   PetscReal      lval, rval, rtmp; /* Terms in the quadature sum to the left and right of 0 */
2447:   PetscInt       d;                /* Digits of precision in the integral */
2448:   mpfr_t         pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;

2450:   PetscFunctionBegin;
2451:   PetscCheck(digits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
2452:   /* Create high precision storage */
2453:   mpfr_inits2(PetscCeilReal(safetyFactor * digits * PetscLogReal(10.) / PetscLogReal(2.)), alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
2454:   /* Initialization */
2455:   mpfr_set_d(alpha, 0.5 * (b - a), MPFR_RNDN);
2456:   mpfr_set_d(beta, 0.5 * (b + a), MPFR_RNDN);
2457:   mpfr_set_d(osum, 0.0, MPFR_RNDN);
2458:   mpfr_set_d(psum, 0.0, MPFR_RNDN);
2459:   mpfr_set_d(h, 1.0, MPFR_RNDN);
2460:   mpfr_const_pi(pi2, MPFR_RNDN);
2461:   mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
2462:   /* Center term */
2463:   rtmp = 0.5 * (b + a);
2464:   func(&rtmp, ctx, &lval);
2465:   mpfr_set(sum, pi2, MPFR_RNDN);
2466:   mpfr_mul(sum, sum, alpha, MPFR_RNDN);
2467:   mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
2468:   /* */
2469:   do {
2470:     PetscReal d1, d2, d3, d4;
2471:     PetscInt  k = 1;

2473:     ++l;
2474:     mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
2475:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %" PetscInt_FMT " sum: %15.15f\n", l, sum); */
2476:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
2477:     mpfr_set(psum, osum, MPFR_RNDN);
2478:     mpfr_set(osum, sum, MPFR_RNDN);
2479:     mpfr_mul_d(h, h, 0.5, MPFR_RNDN);
2480:     mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
2481:     do {
2482:       mpfr_set_si(kh, k, MPFR_RNDN);
2483:       mpfr_mul(kh, kh, h, MPFR_RNDN);
2484:       /* Weight */
2485:       mpfr_set(wk, h, MPFR_RNDN);
2486:       mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
2487:       mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
2488:       mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
2489:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
2490:       mpfr_sqr(tmp, tmp, MPFR_RNDN);
2491:       mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
2492:       mpfr_div(wk, wk, tmp, MPFR_RNDN);
2493:       /* Abscissa */
2494:       mpfr_set_d(yk, 1.0, MPFR_RNDZ);
2495:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
2496:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
2497:       mpfr_exp(tmp, msinh, MPFR_RNDN);
2498:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
2499:       /* Quadrature points */
2500:       mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
2501:       mpfr_mul(lx, lx, alpha, MPFR_RNDU);
2502:       mpfr_add(lx, lx, beta, MPFR_RNDU);
2503:       mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
2504:       mpfr_mul(rx, rx, alpha, MPFR_RNDD);
2505:       mpfr_add(rx, rx, beta, MPFR_RNDD);
2506:       /* Evaluation */
2507:       rtmp = mpfr_get_d(lx, MPFR_RNDU);
2508:       func(&rtmp, ctx, &lval);
2509:       rtmp = mpfr_get_d(rx, MPFR_RNDD);
2510:       func(&rtmp, ctx, &rval);
2511:       /* Update */
2512:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
2513:       mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
2514:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
2515:       mpfr_abs(tmp, tmp, MPFR_RNDN);
2516:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
2517:       mpfr_set(curTerm, tmp, MPFR_RNDN);
2518:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
2519:       mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
2520:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
2521:       mpfr_abs(tmp, tmp, MPFR_RNDN);
2522:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
2523:       mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
2524:       ++k;
2525:       /* Only need to evaluate every other point on refined levels */
2526:       if (l != 1) ++k;
2527:       mpfr_log10(tmp, wk, MPFR_RNDN);
2528:       mpfr_abs(tmp, tmp, MPFR_RNDN);
2529:     } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor * digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
2530:     mpfr_sub(tmp, sum, osum, MPFR_RNDN);
2531:     mpfr_abs(tmp, tmp, MPFR_RNDN);
2532:     mpfr_log10(tmp, tmp, MPFR_RNDN);
2533:     d1 = mpfr_get_d(tmp, MPFR_RNDN);
2534:     mpfr_sub(tmp, sum, psum, MPFR_RNDN);
2535:     mpfr_abs(tmp, tmp, MPFR_RNDN);
2536:     mpfr_log10(tmp, tmp, MPFR_RNDN);
2537:     d2 = mpfr_get_d(tmp, MPFR_RNDN);
2538:     mpfr_log10(tmp, maxTerm, MPFR_RNDN);
2539:     d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
2540:     mpfr_log10(tmp, curTerm, MPFR_RNDN);
2541:     d4 = mpfr_get_d(tmp, MPFR_RNDN);
2542:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1) / d2, 2 * d1), d3), d4)));
2543:   } while (d < digits && l < 8);
2544:   *sol = mpfr_get_d(sum, MPFR_RNDN);
2545:   /* Cleanup */
2546:   mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
2547:   PetscFunctionReturn(PETSC_SUCCESS);
2548: }
2549: #else

2551: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(const PetscReal[], void *, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, void *ctx, PetscReal *sol)
2552: {
2553:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
2554: }
2555: #endif

2557: /*@
2558:   PetscDTTensorQuadratureCreate - create the tensor product quadrature from two lower-dimensional quadratures

2560:   Not Collective

2562:   Input Parameters:
2563: + q1 - The first quadrature
2564: - q2 - The second quadrature

2566:   Output Parameter:
2567: . q - A `PetscQuadrature` object

2569:   Level: intermediate

2571: .seealso: `PetscQuadrature`, `PetscDTGaussTensorQuadrature()`
2572: @*/
2573: PetscErrorCode PetscDTTensorQuadratureCreate(PetscQuadrature q1, PetscQuadrature q2, PetscQuadrature *q)
2574: {
2575:   DMPolytopeType   ct1, ct2, ct;
2576:   const PetscReal *x1, *w1, *x2, *w2;
2577:   PetscReal       *x, *w;
2578:   PetscInt         dim1, Nc1, Np1, order1, qa, d1;
2579:   PetscInt         dim2, Nc2, Np2, order2, qb, d2;
2580:   PetscInt         dim, Nc, Np, order, qc, d;

2582:   PetscFunctionBegin;
2586:   PetscCall(PetscQuadratureGetOrder(q1, &order1));
2587:   PetscCall(PetscQuadratureGetOrder(q2, &order2));
2588:   PetscCheck(order1 == order2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Order1 %" PetscInt_FMT " != %" PetscInt_FMT " Order2", order1, order2);
2589:   PetscCall(PetscQuadratureGetData(q1, &dim1, &Nc1, &Np1, &x1, &w1));
2590:   PetscCall(PetscQuadratureGetCellType(q1, &ct1));
2591:   PetscCall(PetscQuadratureGetData(q2, &dim2, &Nc2, &Np2, &x2, &w2));
2592:   PetscCall(PetscQuadratureGetCellType(q2, &ct2));
2593:   PetscCheck(Nc1 == Nc2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "NumComp1 %" PetscInt_FMT " != %" PetscInt_FMT " NumComp2", Nc1, Nc2);

2595:   switch (ct1) {
2596:   case DM_POLYTOPE_POINT:
2597:     ct = ct2;
2598:     break;
2599:   case DM_POLYTOPE_SEGMENT:
2600:     switch (ct2) {
2601:     case DM_POLYTOPE_POINT:
2602:       ct = ct1;
2603:       break;
2604:     case DM_POLYTOPE_SEGMENT:
2605:       ct = DM_POLYTOPE_QUADRILATERAL;
2606:       break;
2607:     case DM_POLYTOPE_TRIANGLE:
2608:       ct = DM_POLYTOPE_TRI_PRISM;
2609:       break;
2610:     case DM_POLYTOPE_QUADRILATERAL:
2611:       ct = DM_POLYTOPE_HEXAHEDRON;
2612:       break;
2613:     case DM_POLYTOPE_TETRAHEDRON:
2614:       ct = DM_POLYTOPE_UNKNOWN;
2615:       break;
2616:     case DM_POLYTOPE_HEXAHEDRON:
2617:       ct = DM_POLYTOPE_UNKNOWN;
2618:       break;
2619:     default:
2620:       ct = DM_POLYTOPE_UNKNOWN;
2621:     }
2622:     break;
2623:   case DM_POLYTOPE_TRIANGLE:
2624:     switch (ct2) {
2625:     case DM_POLYTOPE_POINT:
2626:       ct = ct1;
2627:       break;
2628:     case DM_POLYTOPE_SEGMENT:
2629:       ct = DM_POLYTOPE_TRI_PRISM;
2630:       break;
2631:     case DM_POLYTOPE_TRIANGLE:
2632:       ct = DM_POLYTOPE_UNKNOWN;
2633:       break;
2634:     case DM_POLYTOPE_QUADRILATERAL:
2635:       ct = DM_POLYTOPE_UNKNOWN;
2636:       break;
2637:     case DM_POLYTOPE_TETRAHEDRON:
2638:       ct = DM_POLYTOPE_UNKNOWN;
2639:       break;
2640:     case DM_POLYTOPE_HEXAHEDRON:
2641:       ct = DM_POLYTOPE_UNKNOWN;
2642:       break;
2643:     default:
2644:       ct = DM_POLYTOPE_UNKNOWN;
2645:     }
2646:     break;
2647:   case DM_POLYTOPE_QUADRILATERAL:
2648:     switch (ct2) {
2649:     case DM_POLYTOPE_POINT:
2650:       ct = ct1;
2651:       break;
2652:     case DM_POLYTOPE_SEGMENT:
2653:       ct = DM_POLYTOPE_HEXAHEDRON;
2654:       break;
2655:     case DM_POLYTOPE_TRIANGLE:
2656:       ct = DM_POLYTOPE_UNKNOWN;
2657:       break;
2658:     case DM_POLYTOPE_QUADRILATERAL:
2659:       ct = DM_POLYTOPE_UNKNOWN;
2660:       break;
2661:     case DM_POLYTOPE_TETRAHEDRON:
2662:       ct = DM_POLYTOPE_UNKNOWN;
2663:       break;
2664:     case DM_POLYTOPE_HEXAHEDRON:
2665:       ct = DM_POLYTOPE_UNKNOWN;
2666:       break;
2667:     default:
2668:       ct = DM_POLYTOPE_UNKNOWN;
2669:     }
2670:     break;
2671:   case DM_POLYTOPE_TETRAHEDRON:
2672:     switch (ct2) {
2673:     case DM_POLYTOPE_POINT:
2674:       ct = ct1;
2675:       break;
2676:     case DM_POLYTOPE_SEGMENT:
2677:       ct = DM_POLYTOPE_UNKNOWN;
2678:       break;
2679:     case DM_POLYTOPE_TRIANGLE:
2680:       ct = DM_POLYTOPE_UNKNOWN;
2681:       break;
2682:     case DM_POLYTOPE_QUADRILATERAL:
2683:       ct = DM_POLYTOPE_UNKNOWN;
2684:       break;
2685:     case DM_POLYTOPE_TETRAHEDRON:
2686:       ct = DM_POLYTOPE_UNKNOWN;
2687:       break;
2688:     case DM_POLYTOPE_HEXAHEDRON:
2689:       ct = DM_POLYTOPE_UNKNOWN;
2690:       break;
2691:     default:
2692:       ct = DM_POLYTOPE_UNKNOWN;
2693:     }
2694:     break;
2695:   case DM_POLYTOPE_HEXAHEDRON:
2696:     switch (ct2) {
2697:     case DM_POLYTOPE_POINT:
2698:       ct = ct1;
2699:       break;
2700:     case DM_POLYTOPE_SEGMENT:
2701:       ct = DM_POLYTOPE_UNKNOWN;
2702:       break;
2703:     case DM_POLYTOPE_TRIANGLE:
2704:       ct = DM_POLYTOPE_UNKNOWN;
2705:       break;
2706:     case DM_POLYTOPE_QUADRILATERAL:
2707:       ct = DM_POLYTOPE_UNKNOWN;
2708:       break;
2709:     case DM_POLYTOPE_TETRAHEDRON:
2710:       ct = DM_POLYTOPE_UNKNOWN;
2711:       break;
2712:     case DM_POLYTOPE_HEXAHEDRON:
2713:       ct = DM_POLYTOPE_UNKNOWN;
2714:       break;
2715:     default:
2716:       ct = DM_POLYTOPE_UNKNOWN;
2717:     }
2718:     break;
2719:   default:
2720:     ct = DM_POLYTOPE_UNKNOWN;
2721:   }
2722:   dim   = dim1 + dim2;
2723:   Nc    = Nc1;
2724:   Np    = Np1 * Np2;
2725:   order = order1;
2726:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, q));
2727:   PetscCall(PetscQuadratureSetCellType(*q, ct));
2728:   PetscCall(PetscQuadratureSetOrder(*q, order));
2729:   PetscCall(PetscMalloc1(Np * dim, &x));
2730:   PetscCall(PetscMalloc1(Np, &w));
2731:   for (qa = 0, qc = 0; qa < Np1; ++qa) {
2732:     for (qb = 0; qb < Np2; ++qb, ++qc) {
2733:       for (d1 = 0, d = 0; d1 < dim1; ++d1, ++d) x[qc * dim + d] = x1[qa * dim1 + d1];
2734:       for (d2 = 0; d2 < dim2; ++d2, ++d) x[qc * dim + d] = x2[qb * dim2 + d2];
2735:       w[qc] = w1[qa] * w2[qb];
2736:     }
2737:   }
2738:   PetscCall(PetscQuadratureSetData(*q, dim, Nc, Np, x, w));
2739:   PetscFunctionReturn(PETSC_SUCCESS);
2740: }

2742: /* Overwrites A. Can only handle full-rank problems with m>=n
2743:    A in column-major format
2744:    Ainv in row-major format
2745:    tau has length m
2746:    worksize must be >= max(1,n)
2747:  */
2748: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m, PetscInt mstride, PetscInt n, PetscReal *A_in, PetscReal *Ainv_out, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2749: {
2750:   PetscBLASInt M, N, K, lda, ldb, ldwork, info;
2751:   PetscScalar *A, *Ainv, *R, *Q, Alpha;

2753:   PetscFunctionBegin;
2754: #if defined(PETSC_USE_COMPLEX)
2755:   {
2756:     PetscInt i, j;
2757:     PetscCall(PetscMalloc2(m * n, &A, m * n, &Ainv));
2758:     for (j = 0; j < n; j++) {
2759:       for (i = 0; i < m; i++) A[i + m * j] = A_in[i + mstride * j];
2760:     }
2761:     mstride = m;
2762:   }
2763: #else
2764:   A    = A_in;
2765:   Ainv = Ainv_out;
2766: #endif

2768:   PetscCall(PetscBLASIntCast(m, &M));
2769:   PetscCall(PetscBLASIntCast(n, &N));
2770:   PetscCall(PetscBLASIntCast(mstride, &lda));
2771:   PetscCall(PetscBLASIntCast(worksize, &ldwork));
2772:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2773:   PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
2774:   PetscCall(PetscFPTrapPop());
2775:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
2776:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

2778:   /* Extract an explicit representation of Q */
2779:   Q = Ainv;
2780:   PetscCall(PetscArraycpy(Q, A, mstride * n));
2781:   K = N; /* full rank */
2782:   PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
2783:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");

2785:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2786:   Alpha = 1.0;
2787:   ldb   = lda;
2788:   PetscCallBLAS("BLAStrsm", BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb));
2789:   /* Ainv is Q, overwritten with inverse */

2791: #if defined(PETSC_USE_COMPLEX)
2792:   {
2793:     PetscInt i;
2794:     for (i = 0; i < m * n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
2795:     PetscCall(PetscFree2(A, Ainv));
2796:   }
2797: #endif
2798:   PetscFunctionReturn(PETSC_SUCCESS);
2799: }

2801: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
2802: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval, const PetscReal *x, PetscInt ndegree, const PetscInt *degrees, PetscBool Transpose, PetscReal *B)
2803: {
2804:   PetscReal *Bv;
2805:   PetscInt   i, j;

2807:   PetscFunctionBegin;
2808:   PetscCall(PetscMalloc1((ninterval + 1) * ndegree, &Bv));
2809:   /* Point evaluation of L_p on all the source vertices */
2810:   PetscCall(PetscDTLegendreEval(ninterval + 1, x, ndegree, degrees, Bv, NULL, NULL));
2811:   /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
2812:   for (i = 0; i < ninterval; i++) {
2813:     for (j = 0; j < ndegree; j++) {
2814:       if (Transpose) B[i + ninterval * j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j];
2815:       else B[i * ndegree + j] = Bv[(i + 1) * ndegree + j] - Bv[i * ndegree + j];
2816:     }
2817:   }
2818:   PetscCall(PetscFree(Bv));
2819:   PetscFunctionReturn(PETSC_SUCCESS);
2820: }

2822: /*@
2823:    PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals

2825:    Not Collective

2827:    Input Parameters:
2828: +  degree - degree of reconstruction polynomial
2829: .  nsource - number of source intervals
2830: .  sourcex - sorted coordinates of source cell boundaries (length nsource+1)
2831: .  ntarget - number of target intervals
2832: -  targetx - sorted coordinates of target cell boundaries (length ntarget+1)

2834:    Output Parameter:
2835: .  R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]

2837:    Level: advanced

2839: .seealso: `PetscDTLegendreEval()`
2840: @*/
2841: PetscErrorCode PetscDTReconstructPoly(PetscInt degree, PetscInt nsource, const PetscReal *sourcex, PetscInt ntarget, const PetscReal *targetx, PetscReal *R)
2842: {
2843:   PetscInt     i, j, k, *bdegrees, worksize;
2844:   PetscReal    xmin, xmax, center, hscale, *sourcey, *targety, *Bsource, *Bsinv, *Btarget;
2845:   PetscScalar *tau, *work;

2847:   PetscFunctionBegin;
2851:   PetscCheck(degree < nsource, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Reconstruction degree %" PetscInt_FMT " must be less than number of source intervals %" PetscInt_FMT, degree, nsource);
2852:   if (PetscDefined(USE_DEBUG)) {
2853:     for (i = 0; i < nsource; i++) PetscCheck(sourcex[i] < sourcex[i + 1], PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Source interval %" PetscInt_FMT " has negative orientation (%g,%g)", i, (double)sourcex[i], (double)sourcex[i + 1]);
2854:     for (i = 0; i < ntarget; i++) PetscCheck(targetx[i] < targetx[i + 1], PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Target interval %" PetscInt_FMT " has negative orientation (%g,%g)", i, (double)targetx[i], (double)targetx[i + 1]);
2855:   }
2856:   xmin     = PetscMin(sourcex[0], targetx[0]);
2857:   xmax     = PetscMax(sourcex[nsource], targetx[ntarget]);
2858:   center   = (xmin + xmax) / 2;
2859:   hscale   = (xmax - xmin) / 2;
2860:   worksize = nsource;
2861:   PetscCall(PetscMalloc4(degree + 1, &bdegrees, nsource + 1, &sourcey, nsource * (degree + 1), &Bsource, worksize, &work));
2862:   PetscCall(PetscMalloc4(nsource, &tau, nsource * (degree + 1), &Bsinv, ntarget + 1, &targety, ntarget * (degree + 1), &Btarget));
2863:   for (i = 0; i <= nsource; i++) sourcey[i] = (sourcex[i] - center) / hscale;
2864:   for (i = 0; i <= degree; i++) bdegrees[i] = i + 1;
2865:   PetscCall(PetscDTLegendreIntegrate(nsource, sourcey, degree + 1, bdegrees, PETSC_TRUE, Bsource));
2866:   PetscCall(PetscDTPseudoInverseQR(nsource, nsource, degree + 1, Bsource, Bsinv, tau, nsource, work));
2867:   for (i = 0; i <= ntarget; i++) targety[i] = (targetx[i] - center) / hscale;
2868:   PetscCall(PetscDTLegendreIntegrate(ntarget, targety, degree + 1, bdegrees, PETSC_FALSE, Btarget));
2869:   for (i = 0; i < ntarget; i++) {
2870:     PetscReal rowsum = 0;
2871:     for (j = 0; j < nsource; j++) {
2872:       PetscReal sum = 0;
2873:       for (k = 0; k < degree + 1; k++) sum += Btarget[i * (degree + 1) + k] * Bsinv[k * nsource + j];
2874:       R[i * nsource + j] = sum;
2875:       rowsum += sum;
2876:     }
2877:     for (j = 0; j < nsource; j++) R[i * nsource + j] /= rowsum; /* normalize each row */
2878:   }
2879:   PetscCall(PetscFree4(bdegrees, sourcey, Bsource, work));
2880:   PetscCall(PetscFree4(tau, Bsinv, targety, Btarget));
2881:   PetscFunctionReturn(PETSC_SUCCESS);
2882: }

2884: /*@C
2885:    PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points

2887:    Not Collective

2889:    Input Parameters:
2890: +  n - the number of GLL nodes
2891: .  nodes - the GLL nodes
2892: .  weights - the GLL weights
2893: -  f - the function values at the nodes

2895:    Output Parameter:
2896: .  in - the value of the integral

2898:    Level: beginner

2900: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`
2901: @*/
2902: PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n, PetscReal *nodes, PetscReal *weights, const PetscReal *f, PetscReal *in)
2903: {
2904:   PetscInt i;

2906:   PetscFunctionBegin;
2907:   *in = 0.;
2908:   for (i = 0; i < n; i++) *in += f[i] * f[i] * weights[i];
2909:   PetscFunctionReturn(PETSC_SUCCESS);
2910: }

2912: /*@C
2913:    PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element

2915:    Not Collective

2917:    Input Parameters:
2918: +  n - the number of GLL nodes
2919: .  nodes - the GLL nodes
2920: -  weights - the GLL weights

2922:    Output Parameter:
2923: .  A - the stiffness element

2925:    Level: beginner

2927:    Notes:
2928:    Destroy this with `PetscGaussLobattoLegendreElementLaplacianDestroy()`

2930:    You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented (the array is symmetric)

2932: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`
2933: @*/
2934: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
2935: {
2936:   PetscReal      **A;
2937:   const PetscReal *gllnodes = nodes;
2938:   const PetscInt   p        = n - 1;
2939:   PetscReal        z0, z1, z2 = -1, x, Lpj, Lpr;
2940:   PetscInt         i, j, nn, r;

2942:   PetscFunctionBegin;
2943:   PetscCall(PetscMalloc1(n, &A));
2944:   PetscCall(PetscMalloc1(n * n, &A[0]));
2945:   for (i = 1; i < n; i++) A[i] = A[i - 1] + n;

2947:   for (j = 1; j < p; j++) {
2948:     x  = gllnodes[j];
2949:     z0 = 1.;
2950:     z1 = x;
2951:     for (nn = 1; nn < p; nn++) {
2952:       z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2953:       z0 = z1;
2954:       z1 = z2;
2955:     }
2956:     Lpj = z2;
2957:     for (r = 1; r < p; r++) {
2958:       if (r == j) {
2959:         A[j][j] = 2. / (3. * (1. - gllnodes[j] * gllnodes[j]) * Lpj * Lpj);
2960:       } else {
2961:         x  = gllnodes[r];
2962:         z0 = 1.;
2963:         z1 = x;
2964:         for (nn = 1; nn < p; nn++) {
2965:           z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2966:           z0 = z1;
2967:           z1 = z2;
2968:         }
2969:         Lpr     = z2;
2970:         A[r][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * Lpr * (gllnodes[j] - gllnodes[r]) * (gllnodes[j] - gllnodes[r]));
2971:       }
2972:     }
2973:   }
2974:   for (j = 1; j < p + 1; j++) {
2975:     x  = gllnodes[j];
2976:     z0 = 1.;
2977:     z1 = x;
2978:     for (nn = 1; nn < p; nn++) {
2979:       z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2980:       z0 = z1;
2981:       z1 = z2;
2982:     }
2983:     Lpj     = z2;
2984:     A[j][0] = 4. * PetscPowRealInt(-1., p) / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. + gllnodes[j]) * (1. + gllnodes[j]));
2985:     A[0][j] = A[j][0];
2986:   }
2987:   for (j = 0; j < p; j++) {
2988:     x  = gllnodes[j];
2989:     z0 = 1.;
2990:     z1 = x;
2991:     for (nn = 1; nn < p; nn++) {
2992:       z2 = x * z1 * (2. * ((PetscReal)nn) + 1.) / (((PetscReal)nn) + 1.) - z0 * (((PetscReal)nn) / (((PetscReal)nn) + 1.));
2993:       z0 = z1;
2994:       z1 = z2;
2995:     }
2996:     Lpj = z2;

2998:     A[p][j] = 4. / (((PetscReal)p) * (((PetscReal)p) + 1.) * Lpj * (1. - gllnodes[j]) * (1. - gllnodes[j]));
2999:     A[j][p] = A[p][j];
3000:   }
3001:   A[0][0] = 0.5 + (((PetscReal)p) * (((PetscReal)p) + 1.) - 2.) / 6.;
3002:   A[p][p] = A[0][0];
3003:   *AA     = A;
3004:   PetscFunctionReturn(PETSC_SUCCESS);
3005: }

3007: /*@C
3008:    PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element created with `PetscGaussLobattoLegendreElementLaplacianCreate()`

3010:    Not Collective

3012:    Input Parameters:
3013: +  n - the number of GLL nodes
3014: .  nodes - the GLL nodes
3015: .  weights - the GLL weightss
3016: -  A - the stiffness element

3018:    Level: beginner

3020: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`
3021: @*/
3022: PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3023: {
3024:   PetscFunctionBegin;
3025:   PetscCall(PetscFree((*AA)[0]));
3026:   PetscCall(PetscFree(*AA));
3027:   *AA = NULL;
3028:   PetscFunctionReturn(PETSC_SUCCESS);
3029: }

3031: /*@C
3032:    PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element

3034:    Not Collective

3036:    Input Parameter:
3037: +  n - the number of GLL nodes
3038: .  nodes - the GLL nodes
3039: .  weights - the GLL weights

3041:    Output Parameters:
3042: .  AA - the stiffness element
3043: -  AAT - the transpose of AA (pass in `NULL` if you do not need this array)

3045:    Level: beginner

3047:    Notes:
3048:    Destroy this with `PetscGaussLobattoLegendreElementGradientDestroy()`

3050:    You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented

3052: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianDestroy()`, `PetscGaussLobattoLegendreElementGradientDestroy()`
3053: @*/
3054: PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA, PetscReal ***AAT)
3055: {
3056:   PetscReal      **A, **AT = NULL;
3057:   const PetscReal *gllnodes = nodes;
3058:   const PetscInt   p        = n - 1;
3059:   PetscReal        Li, Lj, d0;
3060:   PetscInt         i, j;

3062:   PetscFunctionBegin;
3063:   PetscCall(PetscMalloc1(n, &A));
3064:   PetscCall(PetscMalloc1(n * n, &A[0]));
3065:   for (i = 1; i < n; i++) A[i] = A[i - 1] + n;

3067:   if (AAT) {
3068:     PetscCall(PetscMalloc1(n, &AT));
3069:     PetscCall(PetscMalloc1(n * n, &AT[0]));
3070:     for (i = 1; i < n; i++) AT[i] = AT[i - 1] + n;
3071:   }

3073:   if (n == 1) A[0][0] = 0.;
3074:   d0 = (PetscReal)p * ((PetscReal)p + 1.) / 4.;
3075:   for (i = 0; i < n; i++) {
3076:     for (j = 0; j < n; j++) {
3077:       A[i][j] = 0.;
3078:       PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li));
3079:       PetscCall(PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj));
3080:       if (i != j) A[i][j] = Li / (Lj * (gllnodes[i] - gllnodes[j]));
3081:       if ((j == i) && (i == 0)) A[i][j] = -d0;
3082:       if (j == i && i == p) A[i][j] = d0;
3083:       if (AT) AT[j][i] = A[i][j];
3084:     }
3085:   }
3086:   if (AAT) *AAT = AT;
3087:   *AA = A;
3088:   PetscFunctionReturn(PETSC_SUCCESS);
3089: }

3091: /*@C
3092:    PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with `PetscGaussLobattoLegendreElementGradientCreate()`

3094:    Not Collective

3096:    Input Parameters:
3097: +  n - the number of GLL nodes
3098: .  nodes - the GLL nodes
3099: .  weights - the GLL weights
3100: .  AA - the stiffness element
3101: -  AAT - the transpose of the element

3103:    Level: beginner

3105: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionCreate()`
3106: @*/
3107: PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA, PetscReal ***AAT)
3108: {
3109:   PetscFunctionBegin;
3110:   PetscCall(PetscFree((*AA)[0]));
3111:   PetscCall(PetscFree(*AA));
3112:   *AA = NULL;
3113:   if (AAT) {
3114:     PetscCall(PetscFree((*AAT)[0]));
3115:     PetscCall(PetscFree(*AAT));
3116:     *AAT = NULL;
3117:   }
3118:   PetscFunctionReturn(PETSC_SUCCESS);
3119: }

3121: /*@C
3122:    PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element

3124:    Not Collective

3126:    Input Parameters:
3127: +  n - the number of GLL nodes
3128: .  nodes - the GLL nodes
3129: -  weights - the GLL weightss

3131:    Output Parameter:
3132: .  AA - the stiffness element

3134:    Level: beginner

3136:    Notes:
3137:    Destroy this with `PetscGaussLobattoLegendreElementAdvectionDestroy()`

3139:    This is the same as the Gradient operator multiplied by the diagonal mass matrix

3141:    You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented

3143: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementLaplacianCreate()`, `PetscGaussLobattoLegendreElementAdvectionDestroy()`
3144: @*/
3145: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3146: {
3147:   PetscReal      **D;
3148:   const PetscReal *gllweights = weights;
3149:   const PetscInt   glln       = n;
3150:   PetscInt         i, j;

3152:   PetscFunctionBegin;
3153:   PetscCall(PetscGaussLobattoLegendreElementGradientCreate(n, nodes, weights, &D, NULL));
3154:   for (i = 0; i < glln; i++) {
3155:     for (j = 0; j < glln; j++) D[i][j] = gllweights[i] * D[i][j];
3156:   }
3157:   *AA = D;
3158:   PetscFunctionReturn(PETSC_SUCCESS);
3159: }

3161: /*@C
3162:    PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element created with `PetscGaussLobattoLegendreElementAdvectionCreate()`

3164:    Not Collective

3166:    Input Parameters:
3167: +  n - the number of GLL nodes
3168: .  nodes - the GLL nodes
3169: .  weights - the GLL weights
3170: -  A - advection

3172:    Level: beginner

3174: .seealso: `PetscDTGaussLobattoLegendreQuadrature()`, `PetscGaussLobattoLegendreElementAdvectionCreate()`
3175: @*/
3176: PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3177: {
3178:   PetscFunctionBegin;
3179:   PetscCall(PetscFree((*AA)[0]));
3180:   PetscCall(PetscFree(*AA));
3181:   *AA = NULL;
3182:   PetscFunctionReturn(PETSC_SUCCESS);
3183: }

3185: PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3186: {
3187:   PetscReal      **A;
3188:   const PetscReal *gllweights = weights;
3189:   const PetscInt   glln       = n;
3190:   PetscInt         i, j;

3192:   PetscFunctionBegin;
3193:   PetscCall(PetscMalloc1(glln, &A));
3194:   PetscCall(PetscMalloc1(glln * glln, &A[0]));
3195:   for (i = 1; i < glln; i++) A[i] = A[i - 1] + glln;
3196:   if (glln == 1) A[0][0] = 0.;
3197:   for (i = 0; i < glln; i++) {
3198:     for (j = 0; j < glln; j++) {
3199:       A[i][j] = 0.;
3200:       if (j == i) A[i][j] = gllweights[i];
3201:     }
3202:   }
3203:   *AA = A;
3204:   PetscFunctionReturn(PETSC_SUCCESS);
3205: }

3207: PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n, PetscReal *nodes, PetscReal *weights, PetscReal ***AA)
3208: {
3209:   PetscFunctionBegin;
3210:   PetscCall(PetscFree((*AA)[0]));
3211:   PetscCall(PetscFree(*AA));
3212:   *AA = NULL;
3213:   PetscFunctionReturn(PETSC_SUCCESS);
3214: }

3216: /*@
3217:   PetscDTIndexToBary - convert an index into a barycentric coordinate.

3219:   Input Parameters:
3220: + len - the desired length of the barycentric tuple (usually 1 more than the dimension it represents, so a barycentric coordinate in a triangle has length 3)
3221: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
3222: - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum)

3224:   Output Parameter:
3225: . coord - will be filled with the barycentric coordinate

3227:   Level: beginner

3229:   Note:
3230:   The indices map to barycentric coordinates in lexicographic order, where the first index is the
3231:   least significant and the last index is the most significant.

3233: .seealso: `PetscDTBaryToIndex()`
3234: @*/
3235: PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[])
3236: {
3237:   PetscInt c, d, s, total, subtotal, nexttotal;

3239:   PetscFunctionBeginHot;
3240:   PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
3241:   PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative");
3242:   if (!len) {
3243:     if (!sum && !index) PetscFunctionReturn(PETSC_SUCCESS);
3244:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
3245:   }
3246:   for (c = 1, total = 1; c <= len; c++) {
3247:     /* total is the number of ways to have a tuple of length c with sum */
3248:     if (index < total) break;
3249:     total = (total * (sum + c)) / c;
3250:   }
3251:   PetscCheck(c <= len, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range");
3252:   for (d = c; d < len; d++) coord[d] = 0;
3253:   for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) {
3254:     /* subtotal is the number of ways to have a tuple of length c with sum s */
3255:     /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */
3256:     if ((index + subtotal) >= total) {
3257:       coord[--c] = sum - s;
3258:       index -= (total - subtotal);
3259:       sum       = s;
3260:       total     = nexttotal;
3261:       subtotal  = 1;
3262:       nexttotal = 1;
3263:       s         = 0;
3264:     } else {
3265:       subtotal  = (subtotal * (c + s)) / (s + 1);
3266:       nexttotal = (nexttotal * (c - 1 + s)) / (s + 1);
3267:       s++;
3268:     }
3269:   }
3270:   PetscFunctionReturn(PETSC_SUCCESS);
3271: }

3273: /*@
3274:   PetscDTBaryToIndex - convert a barycentric coordinate to an index

3276:   Input Parameters:
3277: + len - the desired length of the barycentric tuple (usually 1 more than the dimension it represents, so a barycentric coordinate in a triangle has length 3)
3278: . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
3279: - coord - a barycentric coordinate with the given length and sum

3281:   Output Parameter:
3282: . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum)

3284:   Level: beginner

3286:   Note:
3287:   The indices map to barycentric coordinates in lexicographic order, where the first index is the
3288:   least significant and the last index is the most significant.

3290: .seealso: `PetscDTIndexToBary`
3291: @*/
3292: PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index)
3293: {
3294:   PetscInt c;
3295:   PetscInt i;
3296:   PetscInt total;

3298:   PetscFunctionBeginHot;
3299:   PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
3300:   if (!len) {
3301:     if (!sum) {
3302:       *index = 0;
3303:       PetscFunctionReturn(PETSC_SUCCESS);
3304:     }
3305:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
3306:   }
3307:   for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c;
3308:   i = total - 1;
3309:   c = len - 1;
3310:   sum -= coord[c];
3311:   while (sum > 0) {
3312:     PetscInt subtotal;
3313:     PetscInt s;

3315:     for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s;
3316:     i -= subtotal;
3317:     sum -= coord[--c];
3318:   }
3319:   *index = i;
3320:   PetscFunctionReturn(PETSC_SUCCESS);
3321: }

3323: /*@
3324:   PetscQuadratureComputePermutations - Compute permutations of quadrature points corresponding to domain orientations

3326:   Input Parameter:
3327: . quad - The `PetscQuadrature`

3329:   Output Parameters:
3330: + Np   - The number of domain orientations
3331: - perm - An array of `IS` permutations, one for ech orientation,

3333:   Level: developer

3335: .seealso: `PetscQuadratureSetCellType()`, `PetscQuadrature`
3336: @*/
3337: PetscErrorCode PetscQuadratureComputePermutations(PetscQuadrature quad, PetscInt *Np, IS *perm[])
3338: {
3339:   DMPolytopeType   ct;
3340:   const PetscReal *xq, *wq;
3341:   PetscInt         dim, qdim, d, Na, o, Nq, q, qp;

3343:   PetscFunctionBegin;
3344:   PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &xq, &wq));
3345:   PetscCall(PetscQuadratureGetCellType(quad, &ct));
3346:   dim = DMPolytopeTypeGetDim(ct);
3347:   Na  = DMPolytopeTypeGetNumArrangments(ct);
3348:   PetscCall(PetscMalloc1(Na, perm));
3349:   if (Np) *Np = Na;
3350:   Na /= 2;
3351:   for (o = -Na; o < Na; ++o) {
3352:     DM        refdm;
3353:     PetscInt *idx;
3354:     PetscReal xi0[3] = {-1., -1., -1.}, v0[3], J[9], detJ, txq[3];
3355:     PetscBool flg;

3357:     PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
3358:     PetscCall(DMPlexOrientPoint(refdm, 0, o));
3359:     PetscCall(DMPlexComputeCellGeometryFEM(refdm, 0, NULL, v0, J, NULL, &detJ));
3360:     PetscCall(PetscMalloc1(Nq, &idx));
3361:     for (q = 0; q < Nq; ++q) {
3362:       CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q * dim], txq);
3363:       for (qp = 0; qp < Nq; ++qp) {
3364:         PetscReal diff = 0.;

3366:         for (d = 0; d < dim; ++d) diff += PetscAbsReal(txq[d] - xq[qp * dim + d]);
3367:         if (diff < PETSC_SMALL) break;
3368:       }
3369:       PetscCheck(qp < Nq, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Transformed quad point %" PetscInt_FMT " does not match another quad point", q);
3370:       idx[q] = qp;
3371:     }
3372:     PetscCall(DMDestroy(&refdm));
3373:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nq, idx, PETSC_OWN_POINTER, &(*perm)[o + Na]));
3374:     PetscCall(ISGetInfo((*perm)[o + Na], IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg));
3375:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ordering for orientation %" PetscInt_FMT " was not a permutation", o);
3376:     PetscCall(ISSetPermutation((*perm)[o + Na]));
3377:   }
3378:   if (!Na) (*perm)[0] = NULL;
3379:   PetscFunctionReturn(PETSC_SUCCESS);
3380: }

3382: /*@
3383:   PetscDTCreateDefaultQuadrature - Create default quadrature for a given cell

3385:   Not collective

3387:   Input Parameters:
3388: + ct     - The integration domain
3389: - qorder - The desired quadrature order

3391:   Output Parameters:
3392: + q  - The cell quadrature
3393: - fq - The face quadrature

3395:   Level: developer

3397: .seealso: `PetscFECreateDefault()`, `PetscDTGaussTensorQuadrature()`, `PetscDTSimplexQuadrature()`, `PetscDTTensorQuadratureCreate()`
3398: @*/
3399: PetscErrorCode PetscDTCreateDefaultQuadrature(DMPolytopeType ct, PetscInt qorder, PetscQuadrature *q, PetscQuadrature *fq)
3400: {
3401:   const PetscInt quadPointsPerEdge = PetscMax(qorder + 1, 1);
3402:   const PetscInt dim               = DMPolytopeTypeGetDim(ct);

3404:   PetscFunctionBegin;
3405:   switch (ct) {
3406:   case DM_POLYTOPE_SEGMENT:
3407:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
3408:   case DM_POLYTOPE_QUADRILATERAL:
3409:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
3410:   case DM_POLYTOPE_HEXAHEDRON:
3411:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
3412:     PetscCall(PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, q));
3413:     PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, fq));
3414:     break;
3415:   case DM_POLYTOPE_TRIANGLE:
3416:   case DM_POLYTOPE_TETRAHEDRON:
3417:     PetscCall(PetscDTSimplexQuadrature(dim, 2 * qorder, PETSCDTSIMPLEXQUAD_DEFAULT, q));
3418:     PetscCall(PetscDTSimplexQuadrature(dim - 1, 2 * qorder, PETSCDTSIMPLEXQUAD_DEFAULT, fq));
3419:     break;
3420:   case DM_POLYTOPE_TRI_PRISM:
3421:   case DM_POLYTOPE_TRI_PRISM_TENSOR: {
3422:     PetscQuadrature q1, q2;

3424:     // TODO: this should be able to use symmetric rules, but doing so causes tests to fail
3425:     PetscCall(PetscDTSimplexQuadrature(2, 2 * qorder, PETSCDTSIMPLEXQUAD_CONIC, &q1));
3426:     PetscCall(PetscDTGaussTensorQuadrature(1, 1, quadPointsPerEdge, -1.0, 1.0, &q2));
3427:     PetscCall(PetscDTTensorQuadratureCreate(q1, q2, q));
3428:     PetscCall(PetscQuadratureDestroy(&q2));
3429:     *fq = q1;
3430:     /* TODO Need separate quadratures for each face */
3431:   } break;
3432:   default:
3433:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No quadrature for celltype %s", DMPolytopeTypes[PetscMin(ct, DM_POLYTOPE_UNKNOWN)]);
3434:   }
3435:   PetscFunctionReturn(PETSC_SUCCESS);
3436: }