1 /* Discretization tools */
2
3 #include <petscdt.h> /*I "petscdt.h" I*/
4 #include <petscblaslapack.h>
5 #include <petsc/private/petscimpl.h>
6 #include <petsc/private/dtimpl.h>
7 #include <petscviewer.h>
8 #include <petscdmplex.h>
9 #include <petscdmshell.h>
10
11 #if defined(PETSC_HAVE_MPFR)
12 #include <mpfr.h>
13 #endif
14
15 const char *const PetscDTNodeTypes[] = {"gaussjacobi", "equispaced", "tanhsinh", "PETSCDTNODES_", NULL};
16
17 static PetscBool GolubWelschCite = PETSC_FALSE;
18 const char GolubWelschCitation[] = "@article{GolubWelsch1969,\n"
19 " author = {Golub and Welsch},\n"
20 " title = {Calculation of Quadrature Rules},\n"
21 " journal = {Math. Comp.},\n"
22 " volume = {23},\n"
23 " number = {106},\n"
24 " pages = {221--230},\n"
25 " year = {1969}\n}\n";
26
27 /* Numerical tests in src/dm/dt/tests/ex1.c show that when computing the nodes and weights of Gauss-Jacobi
28 quadrature rules:
29
30 - in double precision, Newton's method and Golub & Welsch both work for moderate degrees (< 100),
31 - in single precision, Newton's method starts producing incorrect roots around n = 15, but
32 the weights from Golub & Welsch become a problem before then: they produces errors
33 in computing the Jacobi-polynomial Gram matrix around n = 6.
34
35 So we default to Newton's method (required fewer dependencies) */
36 PetscBool PetscDTGaussQuadratureNewton_Internal = PETSC_TRUE;
37
38 PetscClassId PETSCQUADRATURE_CLASSID = 0;
39
40 /*@
41 PetscQuadratureCreate - Create a PetscQuadrature object
42
43 Collective
44
45 Input Parameter:
46 . comm - The communicator for the PetscQuadrature object
47
48 Output Parameter:
49 . q - The PetscQuadrature object
50
51 Level: beginner
52
53 .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
54 @*/
PetscQuadratureCreate(MPI_Comm comm,PetscQuadrature * q)55 PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
56 {
57 PetscErrorCode ierr;
58
59 PetscFunctionBegin;
60 PetscValidPointer(q, 2);
61 ierr = DMInitializePackage();CHKERRQ(ierr);
62 ierr = PetscHeaderCreate(*q,PETSCQUADRATURE_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);CHKERRQ(ierr);
63 (*q)->dim = -1;
64 (*q)->Nc = 1;
65 (*q)->order = -1;
66 (*q)->numPoints = 0;
67 (*q)->points = NULL;
68 (*q)->weights = NULL;
69 PetscFunctionReturn(0);
70 }
71
72 /*@
73 PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object
74
75 Collective on q
76
77 Input Parameter:
78 . q - The PetscQuadrature object
79
80 Output Parameter:
81 . r - The new PetscQuadrature object
82
83 Level: beginner
84
85 .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData()
86 @*/
PetscQuadratureDuplicate(PetscQuadrature q,PetscQuadrature * r)87 PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
88 {
89 PetscInt order, dim, Nc, Nq;
90 const PetscReal *points, *weights;
91 PetscReal *p, *w;
92 PetscErrorCode ierr;
93
94 PetscFunctionBegin;
95 PetscValidPointer(q, 2);
96 ierr = PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);CHKERRQ(ierr);
97 ierr = PetscQuadratureGetOrder(q, &order);CHKERRQ(ierr);
98 ierr = PetscQuadratureSetOrder(*r, order);CHKERRQ(ierr);
99 ierr = PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);CHKERRQ(ierr);
100 ierr = PetscMalloc1(Nq*dim, &p);CHKERRQ(ierr);
101 ierr = PetscMalloc1(Nq*Nc, &w);CHKERRQ(ierr);
102 ierr = PetscArraycpy(p, points, Nq*dim);CHKERRQ(ierr);
103 ierr = PetscArraycpy(w, weights, Nc * Nq);CHKERRQ(ierr);
104 ierr = PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);CHKERRQ(ierr);
105 PetscFunctionReturn(0);
106 }
107
108 /*@
109 PetscQuadratureDestroy - Destroys a PetscQuadrature object
110
111 Collective on q
112
113 Input Parameter:
114 . q - The PetscQuadrature object
115
116 Level: beginner
117
118 .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
119 @*/
PetscQuadratureDestroy(PetscQuadrature * q)120 PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
121 {
122 PetscErrorCode ierr;
123
124 PetscFunctionBegin;
125 if (!*q) PetscFunctionReturn(0);
126 PetscValidHeaderSpecific((*q),PETSCQUADRATURE_CLASSID,1);
127 if (--((PetscObject)(*q))->refct > 0) {
128 *q = NULL;
129 PetscFunctionReturn(0);
130 }
131 ierr = PetscFree((*q)->points);CHKERRQ(ierr);
132 ierr = PetscFree((*q)->weights);CHKERRQ(ierr);
133 ierr = PetscHeaderDestroy(q);CHKERRQ(ierr);
134 PetscFunctionReturn(0);
135 }
136
137 /*@
138 PetscQuadratureGetOrder - Return the order of the method
139
140 Not collective
141
142 Input Parameter:
143 . q - The PetscQuadrature object
144
145 Output Parameter:
146 . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
147
148 Level: intermediate
149
150 .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
151 @*/
PetscQuadratureGetOrder(PetscQuadrature q,PetscInt * order)152 PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
153 {
154 PetscFunctionBegin;
155 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1);
156 PetscValidPointer(order, 2);
157 *order = q->order;
158 PetscFunctionReturn(0);
159 }
160
161 /*@
162 PetscQuadratureSetOrder - Return the order of the method
163
164 Not collective
165
166 Input Parameters:
167 + q - The PetscQuadrature object
168 - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
169
170 Level: intermediate
171
172 .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
173 @*/
PetscQuadratureSetOrder(PetscQuadrature q,PetscInt order)174 PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
175 {
176 PetscFunctionBegin;
177 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1);
178 q->order = order;
179 PetscFunctionReturn(0);
180 }
181
182 /*@
183 PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated
184
185 Not collective
186
187 Input Parameter:
188 . q - The PetscQuadrature object
189
190 Output Parameter:
191 . Nc - The number of components
192
193 Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
194
195 Level: intermediate
196
197 .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
198 @*/
PetscQuadratureGetNumComponents(PetscQuadrature q,PetscInt * Nc)199 PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
200 {
201 PetscFunctionBegin;
202 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1);
203 PetscValidPointer(Nc, 2);
204 *Nc = q->Nc;
205 PetscFunctionReturn(0);
206 }
207
208 /*@
209 PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated
210
211 Not collective
212
213 Input Parameters:
214 + q - The PetscQuadrature object
215 - Nc - The number of components
216
217 Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
218
219 Level: intermediate
220
221 .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
222 @*/
PetscQuadratureSetNumComponents(PetscQuadrature q,PetscInt Nc)223 PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
224 {
225 PetscFunctionBegin;
226 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1);
227 q->Nc = Nc;
228 PetscFunctionReturn(0);
229 }
230
231 /*@C
232 PetscQuadratureGetData - Returns the data defining the quadrature
233
234 Not collective
235
236 Input Parameter:
237 . q - The PetscQuadrature object
238
239 Output Parameters:
240 + dim - The spatial dimension
241 . Nc - The number of components
242 . npoints - The number of quadrature points
243 . points - The coordinates of each quadrature point
244 - weights - The weight of each quadrature point
245
246 Level: intermediate
247
248 Fortran Notes:
249 From Fortran you must call PetscQuadratureRestoreData() when you are done with the data
250
251 .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
252 @*/
PetscQuadratureGetData(PetscQuadrature q,PetscInt * dim,PetscInt * Nc,PetscInt * npoints,const PetscReal * points[],const PetscReal * weights[])253 PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
254 {
255 PetscFunctionBegin;
256 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1);
257 if (dim) {
258 PetscValidPointer(dim, 2);
259 *dim = q->dim;
260 }
261 if (Nc) {
262 PetscValidPointer(Nc, 3);
263 *Nc = q->Nc;
264 }
265 if (npoints) {
266 PetscValidPointer(npoints, 4);
267 *npoints = q->numPoints;
268 }
269 if (points) {
270 PetscValidPointer(points, 5);
271 *points = q->points;
272 }
273 if (weights) {
274 PetscValidPointer(weights, 6);
275 *weights = q->weights;
276 }
277 PetscFunctionReturn(0);
278 }
279
PetscDTJacobianInverse_Internal(PetscInt m,PetscInt n,const PetscReal J[],PetscReal Jinv[])280 static PetscErrorCode PetscDTJacobianInverse_Internal(PetscInt m, PetscInt n, const PetscReal J[], PetscReal Jinv[])
281 {
282 PetscScalar *Js, *Jinvs;
283 PetscInt i, j, k;
284 PetscBLASInt bm, bn, info;
285 PetscErrorCode ierr;
286
287 PetscFunctionBegin;
288 if (!m || !n) PetscFunctionReturn(0);
289 ierr = PetscBLASIntCast(m, &bm);CHKERRQ(ierr);
290 ierr = PetscBLASIntCast(n, &bn);CHKERRQ(ierr);
291 #if defined(PETSC_USE_COMPLEX)
292 ierr = PetscMalloc2(m*n, &Js, m*n, &Jinvs);CHKERRQ(ierr);
293 for (i = 0; i < m*n; i++) Js[i] = J[i];
294 #else
295 Js = (PetscReal *) J;
296 Jinvs = Jinv;
297 #endif
298 if (m == n) {
299 PetscBLASInt *pivots;
300 PetscScalar *W;
301
302 ierr = PetscMalloc2(m, &pivots, m, &W);CHKERRQ(ierr);
303
304 ierr = PetscArraycpy(Jinvs, Js, m * m);CHKERRQ(ierr);
305 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, Jinvs, &bm, pivots, &info));
306 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
307 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, Jinvs, &bm, pivots, W, &bm, &info));
308 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
309 ierr = PetscFree2(pivots, W);CHKERRQ(ierr);
310 } else if (m < n) {
311 PetscScalar *JJT;
312 PetscBLASInt *pivots;
313 PetscScalar *W;
314
315 ierr = PetscMalloc1(m*m, &JJT);CHKERRQ(ierr);
316 ierr = PetscMalloc2(m, &pivots, m, &W);CHKERRQ(ierr);
317 for (i = 0; i < m; i++) {
318 for (j = 0; j < m; j++) {
319 PetscScalar val = 0.;
320
321 for (k = 0; k < n; k++) val += Js[i * n + k] * Js[j * n + k];
322 JJT[i * m + j] = val;
323 }
324 }
325
326 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bm, &bm, JJT, &bm, pivots, &info));
327 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
328 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bm, JJT, &bm, pivots, W, &bm, &info));
329 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
330 for (i = 0; i < n; i++) {
331 for (j = 0; j < m; j++) {
332 PetscScalar val = 0.;
333
334 for (k = 0; k < m; k++) val += Js[k * n + i] * JJT[k * m + j];
335 Jinvs[i * m + j] = val;
336 }
337 }
338 ierr = PetscFree2(pivots, W);CHKERRQ(ierr);
339 ierr = PetscFree(JJT);CHKERRQ(ierr);
340 } else {
341 PetscScalar *JTJ;
342 PetscBLASInt *pivots;
343 PetscScalar *W;
344
345 ierr = PetscMalloc1(n*n, &JTJ);CHKERRQ(ierr);
346 ierr = PetscMalloc2(n, &pivots, n, &W);CHKERRQ(ierr);
347 for (i = 0; i < n; i++) {
348 for (j = 0; j < n; j++) {
349 PetscScalar val = 0.;
350
351 for (k = 0; k < m; k++) val += Js[k * n + i] * Js[k * n + j];
352 JTJ[i * n + j] = val;
353 }
354 }
355
356 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bn, &bn, JTJ, &bn, pivots, &info));
357 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
358 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&bn, JTJ, &bn, pivots, W, &bn, &info));
359 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
360 for (i = 0; i < n; i++) {
361 for (j = 0; j < m; j++) {
362 PetscScalar val = 0.;
363
364 for (k = 0; k < n; k++) val += JTJ[i * n + k] * Js[j * n + k];
365 Jinvs[i * m + j] = val;
366 }
367 }
368 ierr = PetscFree2(pivots, W);CHKERRQ(ierr);
369 ierr = PetscFree(JTJ);CHKERRQ(ierr);
370 }
371 #if defined(PETSC_USE_COMPLEX)
372 for (i = 0; i < m*n; i++) Jinv[i] = PetscRealPart(Jinvs[i]);
373 ierr = PetscFree2(Js, Jinvs);CHKERRQ(ierr);
374 #endif
375 PetscFunctionReturn(0);
376 }
377
378 /*@
379 PetscQuadraturePushForward - Push forward a quadrature functional under an affine transformation.
380
381 Collecive on PetscQuadrature
382
383 Input Arguments:
384 + q - the quadrature functional
385 . imageDim - the dimension of the image of the transformation
386 . origin - a point in the original space
387 . originImage - the image of the origin under the transformation
388 . J - the Jacobian of the image: an [imageDim x dim] matrix in row major order
389 - 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]
390
391 Output Arguments:
392 . 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.
393
394 Note: 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.
395
396 Level: intermediate
397
398 .seealso: PetscDTAltVPullback(), PetscDTAltVPullbackMatrix()
399 @*/
PetscQuadraturePushForward(PetscQuadrature q,PetscInt imageDim,const PetscReal origin[],const PetscReal originImage[],const PetscReal J[],PetscInt formDegree,PetscQuadrature * Jinvstarq)400 PetscErrorCode PetscQuadraturePushForward(PetscQuadrature q, PetscInt imageDim, const PetscReal origin[], const PetscReal originImage[], const PetscReal J[], PetscInt formDegree, PetscQuadrature *Jinvstarq)
401 {
402 PetscInt dim, Nc, imageNc, formSize, Ncopies, imageFormSize, Npoints, pt, i, j, c;
403 const PetscReal *points;
404 const PetscReal *weights;
405 PetscReal *imagePoints, *imageWeights;
406 PetscReal *Jinv;
407 PetscReal *Jinvstar;
408 PetscErrorCode ierr;
409
410 PetscFunctionBegin;
411 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1);
412 if (imageDim < PetscAbsInt(formDegree)) SETERRQ2(PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Cannot represent a %D-form in %D dimensions", PetscAbsInt(formDegree), imageDim);
413 ierr = PetscQuadratureGetData(q, &dim, &Nc, &Npoints, &points, &weights);CHKERRQ(ierr);
414 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &formSize);CHKERRQ(ierr);
415 if (Nc % formSize) SETERRQ2(PetscObjectComm((PetscObject)q), PETSC_ERR_ARG_INCOMP, "Number of components %D is not a multiple of formSize %D\n", Nc, formSize);
416 Ncopies = Nc / formSize;
417 ierr = PetscDTBinomialInt(imageDim, PetscAbsInt(formDegree), &imageFormSize);CHKERRQ(ierr);
418 imageNc = Ncopies * imageFormSize;
419 ierr = PetscMalloc1(Npoints * imageDim, &imagePoints);CHKERRQ(ierr);
420 ierr = PetscMalloc1(Npoints * imageNc, &imageWeights);CHKERRQ(ierr);
421 ierr = PetscMalloc2(imageDim * dim, &Jinv, formSize * imageFormSize, &Jinvstar);CHKERRQ(ierr);
422 ierr = PetscDTJacobianInverse_Internal(imageDim, dim, J, Jinv);CHKERRQ(ierr);
423 ierr = PetscDTAltVPullbackMatrix(imageDim, dim, Jinv, formDegree, Jinvstar);CHKERRQ(ierr);
424 for (pt = 0; pt < Npoints; pt++) {
425 const PetscReal *point = &points[pt * dim];
426 PetscReal *imagePoint = &imagePoints[pt * imageDim];
427
428 for (i = 0; i < imageDim; i++) {
429 PetscReal val = originImage[i];
430
431 for (j = 0; j < dim; j++) val += J[i * dim + j] * (point[j] - origin[j]);
432 imagePoint[i] = val;
433 }
434 for (c = 0; c < Ncopies; c++) {
435 const PetscReal *form = &weights[pt * Nc + c * formSize];
436 PetscReal *imageForm = &imageWeights[pt * imageNc + c * imageFormSize];
437
438 for (i = 0; i < imageFormSize; i++) {
439 PetscReal val = 0.;
440
441 for (j = 0; j < formSize; j++) val += Jinvstar[i * formSize + j] * form[j];
442 imageForm[i] = val;
443 }
444 }
445 }
446 ierr = PetscQuadratureCreate(PetscObjectComm((PetscObject)q), Jinvstarq);CHKERRQ(ierr);
447 ierr = PetscQuadratureSetData(*Jinvstarq, imageDim, imageNc, Npoints, imagePoints, imageWeights);CHKERRQ(ierr);
448 ierr = PetscFree2(Jinv, Jinvstar);CHKERRQ(ierr);
449 PetscFunctionReturn(0);
450 }
451
452 /*@C
453 PetscQuadratureSetData - Sets the data defining the quadrature
454
455 Not collective
456
457 Input Parameters:
458 + q - The PetscQuadrature object
459 . dim - The spatial dimension
460 . Nc - The number of components
461 . npoints - The number of quadrature points
462 . points - The coordinates of each quadrature point
463 - weights - The weight of each quadrature point
464
465 Note: This routine owns the references to points and weights, so they must be allocated using PetscMalloc() and the user should not free them.
466
467 Level: intermediate
468
469 .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
470 @*/
PetscQuadratureSetData(PetscQuadrature q,PetscInt dim,PetscInt Nc,PetscInt npoints,const PetscReal points[],const PetscReal weights[])471 PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
472 {
473 PetscFunctionBegin;
474 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1);
475 if (dim >= 0) q->dim = dim;
476 if (Nc >= 0) q->Nc = Nc;
477 if (npoints >= 0) q->numPoints = npoints;
478 if (points) {
479 PetscValidPointer(points, 4);
480 q->points = points;
481 }
482 if (weights) {
483 PetscValidPointer(weights, 5);
484 q->weights = weights;
485 }
486 PetscFunctionReturn(0);
487 }
488
PetscQuadratureView_Ascii(PetscQuadrature quad,PetscViewer v)489 static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
490 {
491 PetscInt q, d, c;
492 PetscViewerFormat format;
493 PetscErrorCode ierr;
494
495 PetscFunctionBegin;
496 if (quad->Nc > 1) {ierr = PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D) with %D components\n", quad->order, quad->numPoints, quad->dim, quad->Nc);CHKERRQ(ierr);}
497 else {ierr = PetscViewerASCIIPrintf(v, "Quadrature of order %D on %D points (dim %D)\n", quad->order, quad->numPoints, quad->dim);CHKERRQ(ierr);}
498 ierr = PetscViewerGetFormat(v, &format);CHKERRQ(ierr);
499 if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0);
500 for (q = 0; q < quad->numPoints; ++q) {
501 ierr = PetscViewerASCIIPrintf(v, "p%D (", q);CHKERRQ(ierr);
502 ierr = PetscViewerASCIIUseTabs(v, PETSC_FALSE);CHKERRQ(ierr);
503 for (d = 0; d < quad->dim; ++d) {
504 if (d) {ierr = PetscViewerASCIIPrintf(v, ", ");CHKERRQ(ierr);}
505 ierr = PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q*quad->dim+d]);CHKERRQ(ierr);
506 }
507 ierr = PetscViewerASCIIPrintf(v, ") ");CHKERRQ(ierr);
508 if (quad->Nc > 1) {ierr = PetscViewerASCIIPrintf(v, "w%D (", q);CHKERRQ(ierr);}
509 for (c = 0; c < quad->Nc; ++c) {
510 if (c) {ierr = PetscViewerASCIIPrintf(v, ", ");CHKERRQ(ierr);}
511 ierr = PetscViewerASCIIPrintf(v, "%+g", (double)quad->weights[q*quad->Nc+c]);CHKERRQ(ierr);
512 }
513 if (quad->Nc > 1) {ierr = PetscViewerASCIIPrintf(v, ")");CHKERRQ(ierr);}
514 ierr = PetscViewerASCIIPrintf(v, "\n");CHKERRQ(ierr);
515 ierr = PetscViewerASCIIUseTabs(v, PETSC_TRUE);CHKERRQ(ierr);
516 }
517 PetscFunctionReturn(0);
518 }
519
520 /*@C
521 PetscQuadratureView - Views a PetscQuadrature object
522
523 Collective on quad
524
525 Input Parameters:
526 + quad - The PetscQuadrature object
527 - viewer - The PetscViewer object
528
529 Level: beginner
530
531 .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
532 @*/
PetscQuadratureView(PetscQuadrature quad,PetscViewer viewer)533 PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
534 {
535 PetscBool iascii;
536 PetscErrorCode ierr;
537
538 PetscFunctionBegin;
539 PetscValidHeader(quad, 1);
540 if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
541 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) quad), &viewer);CHKERRQ(ierr);}
542 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
543 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
544 if (iascii) {ierr = PetscQuadratureView_Ascii(quad, viewer);CHKERRQ(ierr);}
545 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
546 PetscFunctionReturn(0);
547 }
548
549 /*@C
550 PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement
551
552 Not collective
553
554 Input Parameter:
555 + q - The original PetscQuadrature
556 . numSubelements - The number of subelements the original element is divided into
557 . v0 - An array of the initial points for each subelement
558 - jac - An array of the Jacobian mappings from the reference to each subelement
559
560 Output Parameters:
561 . dim - The dimension
562
563 Note: Together v0 and jac define an affine mapping from the original reference element to each subelement
564
565 Not available from Fortran
566
567 Level: intermediate
568
569 .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
570 @*/
PetscQuadratureExpandComposite(PetscQuadrature q,PetscInt numSubelements,const PetscReal v0[],const PetscReal jac[],PetscQuadrature * qref)571 PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
572 {
573 const PetscReal *points, *weights;
574 PetscReal *pointsRef, *weightsRef;
575 PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
576 PetscErrorCode ierr;
577
578 PetscFunctionBegin;
579 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 1);
580 PetscValidPointer(v0, 3);
581 PetscValidPointer(jac, 4);
582 PetscValidPointer(qref, 5);
583 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, qref);CHKERRQ(ierr);
584 ierr = PetscQuadratureGetOrder(q, &order);CHKERRQ(ierr);
585 ierr = PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);CHKERRQ(ierr);
586 npointsRef = npoints*numSubelements;
587 ierr = PetscMalloc1(npointsRef*dim,&pointsRef);CHKERRQ(ierr);
588 ierr = PetscMalloc1(npointsRef*Nc, &weightsRef);CHKERRQ(ierr);
589 for (c = 0; c < numSubelements; ++c) {
590 for (p = 0; p < npoints; ++p) {
591 for (d = 0; d < dim; ++d) {
592 pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
593 for (e = 0; e < dim; ++e) {
594 pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
595 }
596 }
597 /* Could also use detJ here */
598 for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
599 }
600 }
601 ierr = PetscQuadratureSetOrder(*qref, order);CHKERRQ(ierr);
602 ierr = PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);CHKERRQ(ierr);
603 PetscFunctionReturn(0);
604 }
605
606 /* Compute the coefficients for the Jacobi polynomial recurrence,
607 *
608 * J^{a,b}_n(x) = (cnm1 + cnm1x * x) * J^{a,b}_{n-1}(x) - cnm2 * J^{a,b}_{n-2}(x).
609 */
610 #define PetscDTJacobiRecurrence_Internal(n,a,b,cnm1,cnm1x,cnm2) \
611 do { \
612 PetscReal _a = (a); \
613 PetscReal _b = (b); \
614 PetscReal _n = (n); \
615 if (n == 1) { \
616 (cnm1) = (_a-_b) * 0.5; \
617 (cnm1x) = (_a+_b+2.)*0.5; \
618 (cnm2) = 0.; \
619 } else { \
620 PetscReal _2n = _n+_n; \
621 PetscReal _d = (_2n*(_n+_a+_b)*(_2n+_a+_b-2)); \
622 PetscReal _n1 = (_2n+_a+_b-1.)*(_a*_a-_b*_b); \
623 PetscReal _n1x = (_2n+_a+_b-1.)*(_2n+_a+_b)*(_2n+_a+_b-2); \
624 PetscReal _n2 = 2.*((_n+_a-1.)*(_n+_b-1.)*(_2n+_a+_b)); \
625 (cnm1) = _n1 / _d; \
626 (cnm1x) = _n1x / _d; \
627 (cnm2) = _n2 / _d; \
628 } \
629 } while (0)
630
631 /*@
632 PetscDTJacobiNorm - Compute the weighted L2 norm of a Jacobi polynomial.
633
634 $\| P^{\alpha,\beta}_n \|_{\alpha,\beta}^2 = \int_{-1}^1 (1 + x)^{\alpha} (1 - x)^{\beta} P^{\alpha,\beta}_n (x)^2 dx.$
635
636 Input Arguments:
637 - alpha - the left exponent > -1
638 . beta - the right exponent > -1
639 + n - the polynomial degree
640
641 Output Arguments:
642 . norm - the weighted L2 norm
643
644 Level: beginner
645
646 .seealso: PetscDTJacobiEval()
647 @*/
PetscDTJacobiNorm(PetscReal alpha,PetscReal beta,PetscInt n,PetscReal * norm)648 PetscErrorCode PetscDTJacobiNorm(PetscReal alpha, PetscReal beta, PetscInt n, PetscReal *norm)
649 {
650 PetscReal twoab1;
651 PetscReal gr;
652
653 PetscFunctionBegin;
654 if (alpha <= -1.) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent alpha %g <= -1. invalid\n", (double) alpha);
655 if (beta <= -1.) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exponent beta %g <= -1. invalid\n", (double) beta);
656 if (n < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "n %D < 0 invalid\n", n);
657 twoab1 = PetscPowReal(2., alpha + beta + 1.);
658 #if defined(PETSC_HAVE_LGAMMA)
659 if (!n) {
660 gr = PetscExpReal(PetscLGamma(alpha+1.) + PetscLGamma(beta+1.) - PetscLGamma(alpha+beta+2.));
661 } else {
662 gr = PetscExpReal(PetscLGamma(n+alpha+1.) + PetscLGamma(n+beta+1.) - (PetscLGamma(n+1.) + PetscLGamma(n+alpha+beta+1.))) / (n+n+alpha+beta+1.);
663 }
664 #else
665 {
666 PetscInt alphai = (PetscInt) alpha;
667 PetscInt betai = (PetscInt) beta;
668 PetscInt i;
669
670 gr = n ? (1. / (n+n+alpha+beta+1.)) : 1.;
671 if ((PetscReal) alphai == alpha) {
672 if (!n) {
673 for (i = 0; i < alphai; i++) gr *= (i+1.) / (beta+i+1.);
674 gr /= (alpha+beta+1.);
675 } else {
676 for (i = 0; i < alphai; i++) gr *= (n+i+1.) / (n+beta+i+1.);
677 }
678 } else if ((PetscReal) betai == beta) {
679 if (!n) {
680 for (i = 0; i < betai; i++) gr *= (i+1.) / (alpha+i+2.);
681 gr /= (alpha+beta+1.);
682 } else {
683 for (i = 0; i < betai; i++) gr *= (n+i+1.) / (n+alpha+i+1.);
684 }
685 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
686 }
687 #endif
688 *norm = PetscSqrtReal(twoab1 * gr);
689 PetscFunctionReturn(0);
690 }
691
PetscDTJacobiEval_Internal(PetscInt npoints,PetscReal a,PetscReal b,PetscInt k,const PetscReal * points,PetscInt ndegree,const PetscInt * degrees,PetscReal * p)692 static PetscErrorCode PetscDTJacobiEval_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscInt k, const PetscReal *points, PetscInt ndegree, const PetscInt *degrees, PetscReal *p)
693 {
694 PetscReal ak, bk;
695 PetscReal abk1;
696 PetscInt i,l,maxdegree;
697
698 PetscFunctionBegin;
699 maxdegree = degrees[ndegree-1] - k;
700 ak = a + k;
701 bk = b + k;
702 abk1 = a + b + k + 1.;
703 if (maxdegree < 0) {
704 for (i = 0; i < npoints; i++) for (l = 0; l < ndegree; l++) p[i*ndegree+l] = 0.;
705 PetscFunctionReturn(0);
706 }
707 for (i=0; i<npoints; i++) {
708 PetscReal pm1,pm2,x;
709 PetscReal cnm1, cnm1x, cnm2;
710 PetscInt j,m;
711
712 x = points[i];
713 pm2 = 1.;
714 PetscDTJacobiRecurrence_Internal(1,ak,bk,cnm1,cnm1x,cnm2);
715 pm1 = (cnm1 + cnm1x*x);
716 l = 0;
717 while (l < ndegree && degrees[l] - k < 0) {
718 p[l++] = 0.;
719 }
720 while (l < ndegree && degrees[l] - k == 0) {
721 p[l] = pm2;
722 for (m = 0; m < k; m++) p[l] *= (abk1 + m) * 0.5;
723 l++;
724 }
725 while (l < ndegree && degrees[l] - k == 1) {
726 p[l] = pm1;
727 for (m = 0; m < k; m++) p[l] *= (abk1 + 1 + m) * 0.5;
728 l++;
729 }
730 for (j=2; j<=maxdegree; j++) {
731 PetscReal pp;
732
733 PetscDTJacobiRecurrence_Internal(j,ak,bk,cnm1,cnm1x,cnm2);
734 pp = (cnm1 + cnm1x*x)*pm1 - cnm2*pm2;
735 pm2 = pm1;
736 pm1 = pp;
737 while (l < ndegree && degrees[l] - k == j) {
738 p[l] = pp;
739 for (m = 0; m < k; m++) p[l] *= (abk1 + j + m) * 0.5;
740 l++;
741 }
742 }
743 p += ndegree;
744 }
745 PetscFunctionReturn(0);
746 }
747
748 /*@
749 PetscDTJacobiEvalJet - Evaluate the jet (function and derivatives) of the Jacobi polynomials basis up to a given degree. The Jacobi polynomials with indices $\alpha$ and $\beta$ are orthogonal with respect to the weighted inner product $\langle f, g \rangle = \int_{-1}^1 (1+x)^{\alpha} (1-x)^{\beta) f(x) g(x) dx$.
750
751 Input Arguments:
752 + alpha - the left exponent of the weight
753 . beta - the right exponetn of the weight
754 . npoints - the number of points to evaluate the polynomials at
755 . points - [npoints] array of point coordinates
756 . degree - the maximm degree polynomial space to evaluate, (degree + 1) will be evaluated total.
757 - k - the maximum derivative to evaluate in the jet, (k + 1) will be evaluated total.
758
759 Output Argments:
760 - p - an array containing the evaluations of the Jacobi polynomials's jets on the points. the size is (degree + 1) x
761 (k + 1) x npoints, which also describes the order of the dimensions of this three-dimensional array: the first
762 (slowest varying) dimension is polynomial degree; the second dimension is derivative order; the third (fastest
763 varying) dimension is the index of the evaluation point.
764
765 Level: advanced
766
767 .seealso: PetscDTJacobiEval(), PetscDTPKDEvalJet()
768 @*/
PetscDTJacobiEvalJet(PetscReal alpha,PetscReal beta,PetscInt npoints,const PetscReal points[],PetscInt degree,PetscInt k,PetscReal p[])769 PetscErrorCode PetscDTJacobiEvalJet(PetscReal alpha, PetscReal beta, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
770 {
771 PetscInt i, j, l;
772 PetscInt *degrees;
773 PetscReal *psingle;
774 PetscErrorCode ierr;
775
776 PetscFunctionBegin;
777 if (degree == 0) {
778 PetscInt zero = 0;
779
780 for (i = 0; i <= k; i++) {
781 ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, 1, &zero, &p[i*npoints]);CHKERRQ(ierr);
782 }
783 PetscFunctionReturn(0);
784 }
785 ierr = PetscMalloc1(degree + 1, °rees);CHKERRQ(ierr);
786 ierr = PetscMalloc1((degree + 1) * npoints, &psingle);CHKERRQ(ierr);
787 for (i = 0; i <= degree; i++) degrees[i] = i;
788 for (i = 0; i <= k; i++) {
789 ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, i, points, degree + 1, degrees, psingle);CHKERRQ(ierr);
790 for (j = 0; j <= degree; j++) {
791 for (l = 0; l < npoints; l++) {
792 p[(j * (k + 1) + i) * npoints + l] = psingle[l * (degree + 1) + j];
793 }
794 }
795 }
796 ierr = PetscFree(psingle);CHKERRQ(ierr);
797 ierr = PetscFree(degrees);CHKERRQ(ierr);
798 PetscFunctionReturn(0);
799 }
800
801 /*@
802 PetscDTJacobiEval - evaluate Jacobi polynomials for the weight function $(1.+x)^{\alpha} (1.-x)^{\beta}$
803 at points
804
805 Not Collective
806
807 Input Arguments:
808 + npoints - number of spatial points to evaluate at
809 . alpha - the left exponent > -1
810 . beta - the right exponent > -1
811 . points - array of locations to evaluate at
812 . ndegree - number of basis degrees to evaluate
813 - degrees - sorted array of degrees to evaluate
814
815 Output Arguments:
816 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
817 . D - row-oriented derivative evaluation matrix (or NULL)
818 - D2 - row-oriented second derivative evaluation matrix (or NULL)
819
820 Level: intermediate
821
822 .seealso: PetscDTGaussQuadrature()
823 @*/
PetscDTJacobiEval(PetscInt npoints,PetscReal alpha,PetscReal beta,const PetscReal * points,PetscInt ndegree,const PetscInt * degrees,PetscReal * B,PetscReal * D,PetscReal * D2)824 PetscErrorCode PetscDTJacobiEval(PetscInt npoints,PetscReal alpha, PetscReal beta, const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
825 {
826 PetscErrorCode ierr;
827
828 PetscFunctionBegin;
829 if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1.");
830 if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1.");
831 if (!npoints || !ndegree) PetscFunctionReturn(0);
832 if (B) {ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, 0, points, ndegree, degrees, B);CHKERRQ(ierr);}
833 if (D) {ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, 1, points, ndegree, degrees, D);CHKERRQ(ierr);}
834 if (D2) {ierr = PetscDTJacobiEval_Internal(npoints, alpha, beta, 2, points, ndegree, degrees, D2);CHKERRQ(ierr);}
835 PetscFunctionReturn(0);
836 }
837
838 /*@
839 PetscDTLegendreEval - evaluate Legendre polynomials at points
840
841 Not Collective
842
843 Input Arguments:
844 + npoints - number of spatial points to evaluate at
845 . points - array of locations to evaluate at
846 . ndegree - number of basis degrees to evaluate
847 - degrees - sorted array of degrees to evaluate
848
849 Output Arguments:
850 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
851 . D - row-oriented derivative evaluation matrix (or NULL)
852 - D2 - row-oriented second derivative evaluation matrix (or NULL)
853
854 Level: intermediate
855
856 .seealso: PetscDTGaussQuadrature()
857 @*/
PetscDTLegendreEval(PetscInt npoints,const PetscReal * points,PetscInt ndegree,const PetscInt * degrees,PetscReal * B,PetscReal * D,PetscReal * D2)858 PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
859 {
860 PetscErrorCode ierr;
861
862 PetscFunctionBegin;
863 ierr = PetscDTJacobiEval(npoints, 0., 0., points, ndegree, degrees, B, D, D2);CHKERRQ(ierr);
864 PetscFunctionReturn(0);
865 }
866
867 /*@
868 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)
869
870 Input Parameters:
871 + len - the desired length of the degree tuple
872 - index - the index to convert: should be >= 0
873
874 Output Parameter:
875 . degtup - will be filled with a tuple of degrees
876
877 Level: beginner
878
879 Note: for two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
880 acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
881 last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).
882
883 .seealso: PetscDTGradedOrderToIndex()
884 @*/
PetscDTIndexToGradedOrder(PetscInt len,PetscInt index,PetscInt degtup[])885 PetscErrorCode PetscDTIndexToGradedOrder(PetscInt len, PetscInt index, PetscInt degtup[])
886 {
887 PetscInt i, total;
888 PetscInt sum;
889
890 PetscFunctionBeginHot;
891 if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
892 if (index < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative");
893 total = 1;
894 sum = 0;
895 while (index >= total) {
896 index -= total;
897 total = (total * (len + sum)) / (sum + 1);
898 sum++;
899 }
900 for (i = 0; i < len; i++) {
901 PetscInt c;
902
903 degtup[i] = sum;
904 for (c = 0, total = 1; c < sum; c++) {
905 /* going into the loop, total is the number of way to have a tuple of sum exactly c with length len - 1 - i */
906 if (index < total) break;
907 index -= total;
908 total = (total * (len - 1 - i + c)) / (c + 1);
909 degtup[i]--;
910 }
911 sum -= degtup[i];
912 }
913 PetscFunctionReturn(0);
914 }
915
916 /*@
917 PetscDTGradedOrderToIndex - convert a tuple into an index in a graded order, the inverse of PetscDTIndexToGradedOrder().
918
919 Input Parameters:
920 + len - the length of the degree tuple
921 - degtup - tuple with this length
922
923 Output Parameter:
924 . index - index in graded order: >= 0
925
926 Level: Beginner
927
928 Note: for two tuples x and y with the same degree sum, partial degree sums over the final elements of the tuples
929 acts as a tiebreaker. For example, (2, 1, 1) and (1, 2, 1) have the same degree sum, but the degree sum over the
930 last two elements is smaller for the former, so (2, 1, 1) < (1, 2, 1).
931
932 .seealso: PetscDTIndexToGradedOrder()
933 @*/
PetscDTGradedOrderToIndex(PetscInt len,const PetscInt degtup[],PetscInt * index)934 PetscErrorCode PetscDTGradedOrderToIndex(PetscInt len, const PetscInt degtup[], PetscInt *index)
935 {
936 PetscInt i, idx, sum, total;
937
938 PetscFunctionBeginHot;
939 if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
940 for (i = 0, sum = 0; i < len; i++) sum += degtup[i];
941 idx = 0;
942 total = 1;
943 for (i = 0; i < sum; i++) {
944 idx += total;
945 total = (total * (len + i)) / (i + 1);
946 }
947 for (i = 0; i < len - 1; i++) {
948 PetscInt c;
949
950 total = 1;
951 sum -= degtup[i];
952 for (c = 0; c < sum; c++) {
953 idx += total;
954 total = (total * (len - 1 - i + c)) / (c + 1);
955 }
956 }
957 *index = idx;
958 PetscFunctionReturn(0);
959 }
960
961 static PetscBool PKDCite = PETSC_FALSE;
962 const char PKDCitation[] = "@article{Kirby2010,\n"
963 " title={Singularity-free evaluation of collapsed-coordinate orthogonal polynomials},\n"
964 " author={Kirby, Robert C},\n"
965 " journal={ACM Transactions on Mathematical Software (TOMS)},\n"
966 " volume={37},\n"
967 " number={1},\n"
968 " pages={1--16},\n"
969 " year={2010},\n"
970 " publisher={ACM New York, NY, USA}\n}\n";
971
972 /*@
973 PetscDTPKDEvalJet - Evaluate the jet (function and derivatives) of the Prioriol-Koornwinder-Dubiner (PKD) basis for
974 the space of polynomials up to a given degree. The PKD basis is L2-orthonormal on the biunit simplex (which is used
975 as the reference element for finite elements in PETSc), which makes it a stable basis to use for evaluating
976 polynomials in that domain.
977
978 Input Arguments:
979 + dim - the number of variables in the multivariate polynomials
980 . npoints - the number of points to evaluate the polynomials at
981 . points - [npoints x dim] array of point coordinates
982 . 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.
983 - k - the maximum order partial derivative to evaluate in the jet. There are (dim + k choose dim) partial derivatives
984 in the jet. Choosing k = 0 means to evaluate just the function and no derivatives
985
986 Output Argments:
987 - p - an array containing the evaluations of the PKD polynomials' jets on the points. The size is ((dim + degree)
988 choose dim) x ((dim + k) choose dim) x npoints, which also describes the order of the dimensions of this
989 three-dimensional array: the first (slowest varying) dimension is basis function index; the second dimension is jet
990 index; the third (fastest varying) dimension is the index of the evaluation point.
991
992 Level: advanced
993
994 Note: The ordering of the basis functions, and the ordering of the derivatives in the jet, both follow the graded
995 ordering of PetscDTIndexToGradedOrder() and PetscDTGradedOrderToIndex(). For example, in 3D, the polynomial with
996 leading monomial x^3,y^1,z^2, which as degree tuple (2,0,1), which by PetscDTGradedOrderToIndex() has index 12 (it is the 13th basis function in the space);
997 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).
998
999 The implementation uses Kirby's singularity-free evaluation algorithm, https://doi.org/10.1145/1644001.1644006.
1000
1001 .seealso: PetscDTGradedOrderToIndex(), PetscDTIndexToGradedOrder(), PetscDTJacobiEvalJet()
1002 @*/
PetscDTPKDEvalJet(PetscInt dim,PetscInt npoints,const PetscReal points[],PetscInt degree,PetscInt k,PetscReal p[])1003 PetscErrorCode PetscDTPKDEvalJet(PetscInt dim, PetscInt npoints, const PetscReal points[], PetscInt degree, PetscInt k, PetscReal p[])
1004 {
1005 PetscInt degidx, kidx, d, pt;
1006 PetscInt Nk, Ndeg;
1007 PetscInt *ktup, *degtup;
1008 PetscReal *scales, initscale, scaleexp;
1009 PetscErrorCode ierr;
1010
1011 PetscFunctionBegin;
1012 ierr = PetscCitationsRegister(PKDCitation, &PKDCite);CHKERRQ(ierr);
1013 ierr = PetscDTBinomialInt(dim + k, k, &Nk);CHKERRQ(ierr);
1014 ierr = PetscDTBinomialInt(degree + dim, degree, &Ndeg);CHKERRQ(ierr);
1015 ierr = PetscMalloc2(dim, °tup, dim, &ktup);CHKERRQ(ierr);
1016 ierr = PetscMalloc1(Ndeg, &scales);CHKERRQ(ierr);
1017 initscale = 1.;
1018 if (dim > 1) {
1019 ierr = PetscDTBinomial(dim,2,&scaleexp);CHKERRQ(ierr);
1020 initscale = PetscPowReal(2.,scaleexp*0.5);CHKERRQ(ierr);
1021 }
1022 for (degidx = 0; degidx < Ndeg; degidx++) {
1023 PetscInt e, i;
1024 PetscInt m1idx = -1, m2idx = -1;
1025 PetscInt n;
1026 PetscInt degsum;
1027 PetscReal alpha;
1028 PetscReal cnm1, cnm1x, cnm2;
1029 PetscReal norm;
1030
1031 ierr = PetscDTIndexToGradedOrder(dim, degidx, degtup);CHKERRQ(ierr);
1032 for (d = dim - 1; d >= 0; d--) if (degtup[d]) break;
1033 if (d < 0) { /* constant is 1 everywhere, all derivatives are zero */
1034 scales[degidx] = initscale;
1035 for (e = 0; e < dim; e++) {
1036 ierr = PetscDTJacobiNorm(e,0.,0,&norm);CHKERRQ(ierr);
1037 scales[degidx] /= norm;
1038 }
1039 for (i = 0; i < npoints; i++) p[degidx * Nk * npoints + i] = 1.;
1040 for (i = 0; i < (Nk - 1) * npoints; i++) p[(degidx * Nk + 1) * npoints + i] = 0.;
1041 continue;
1042 }
1043 n = degtup[d];
1044 degtup[d]--;
1045 ierr = PetscDTGradedOrderToIndex(dim, degtup, &m1idx);CHKERRQ(ierr);
1046 if (degtup[d] > 0) {
1047 degtup[d]--;
1048 ierr = PetscDTGradedOrderToIndex(dim, degtup, &m2idx);CHKERRQ(ierr);
1049 degtup[d]++;
1050 }
1051 degtup[d]++;
1052 for (e = 0, degsum = 0; e < d; e++) degsum += degtup[e];
1053 alpha = 2 * degsum + d;
1054 PetscDTJacobiRecurrence_Internal(n,alpha,0.,cnm1,cnm1x,cnm2);
1055
1056
1057 scales[degidx] = initscale;
1058 for (e = 0, degsum = 0; e < dim; e++) {
1059 PetscInt f;
1060 PetscReal ealpha;
1061 PetscReal enorm;
1062
1063 ealpha = 2 * degsum + e;
1064 for (f = 0; f < degsum; f++) scales[degidx] *= 2.;
1065 ierr = PetscDTJacobiNorm(ealpha,0.,degtup[e],&enorm);CHKERRQ(ierr);
1066 scales[degidx] /= enorm;
1067 degsum += degtup[e];
1068 }
1069
1070 for (pt = 0; pt < npoints; pt++) {
1071 /* compute the multipliers */
1072 PetscReal thetanm1, thetanm1x, thetanm2;
1073
1074 thetanm1x = dim - (d+1) + 2.*points[pt * dim + d];
1075 for (e = d+1; e < dim; e++) thetanm1x += points[pt * dim + e];
1076 thetanm1x *= 0.5;
1077 thetanm1 = (2. - (dim-(d+1)));
1078 for (e = d+1; e < dim; e++) thetanm1 -= points[pt * dim + e];
1079 thetanm1 *= 0.5;
1080 thetanm2 = thetanm1 * thetanm1;
1081
1082 for (kidx = 0; kidx < Nk; kidx++) {
1083 PetscInt f;
1084
1085 ierr = PetscDTIndexToGradedOrder(dim, kidx, ktup);CHKERRQ(ierr);
1086 /* first sum in the same derivative terms */
1087 p[(degidx * Nk + kidx) * npoints + pt] = (cnm1 * thetanm1 + cnm1x * thetanm1x) * p[(m1idx * Nk + kidx) * npoints + pt];
1088 if (m2idx >= 0) {
1089 p[(degidx * Nk + kidx) * npoints + pt] -= cnm2 * thetanm2 * p[(m2idx * Nk + kidx) * npoints + pt];
1090 }
1091
1092 for (f = d; f < dim; f++) {
1093 PetscInt km1idx, mplty = ktup[f];
1094
1095 if (!mplty) continue;
1096 ktup[f]--;
1097 ierr = PetscDTGradedOrderToIndex(dim, ktup, &km1idx);CHKERRQ(ierr);
1098
1099 /* the derivative of cnm1x * thetanm1x wrt x variable f is 0.5 * cnm1x if f > d otherwise it is cnm1x */
1100 /* the derivative of cnm1 * thetanm1 wrt x variable f is 0 if f == d, otherwise it is -0.5 * cnm1 */
1101 /* the derivative of -cnm2 * thetanm2 wrt x variable f is 0 if f == d, otherwise it is cnm2 * thetanm1 */
1102 if (f > d) {
1103 PetscInt f2;
1104
1105 p[(degidx * Nk + kidx) * npoints + pt] += mplty * 0.5 * (cnm1x - cnm1) * p[(m1idx * Nk + km1idx) * npoints + pt];
1106 if (m2idx >= 0) {
1107 p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm2 * thetanm1 * p[(m2idx * Nk + km1idx) * npoints + pt];
1108 /* second derivatives of -cnm2 * thetanm2 wrt x variable f,f2 is like - 0.5 * cnm2 */
1109 for (f2 = f; f2 < dim; f2++) {
1110 PetscInt km2idx, mplty2 = ktup[f2];
1111 PetscInt factor;
1112
1113 if (!mplty2) continue;
1114 ktup[f2]--;
1115 ierr = PetscDTGradedOrderToIndex(dim, ktup, &km2idx);CHKERRQ(ierr);
1116
1117 factor = mplty * mplty2;
1118 if (f == f2) factor /= 2;
1119 p[(degidx * Nk + kidx) * npoints + pt] -= 0.5 * factor * cnm2 * p[(m2idx * Nk + km2idx) * npoints + pt];
1120 ktup[f2]++;
1121 }
1122 }
1123 } else {
1124 p[(degidx * Nk + kidx) * npoints + pt] += mplty * cnm1x * p[(m1idx * Nk + km1idx) * npoints + pt];
1125 }
1126 ktup[f]++;
1127 }
1128 }
1129 }
1130 }
1131 for (degidx = 0; degidx < Ndeg; degidx++) {
1132 PetscReal scale = scales[degidx];
1133 PetscInt i;
1134
1135 for (i = 0; i < Nk * npoints; i++) p[degidx*Nk*npoints + i] *= scale;
1136 }
1137 ierr = PetscFree(scales);CHKERRQ(ierr);
1138 ierr = PetscFree2(degtup, ktup);CHKERRQ(ierr);
1139 PetscFunctionReturn(0);
1140 }
1141
1142 /* solve the symmetric tridiagonal eigenvalue system, writing the eigenvalues into eigs and the eigenvectors into V
1143 * with lds n; diag and subdiag are overwritten */
PetscDTSymmetricTridiagonalEigensolve(PetscInt n,PetscReal diag[],PetscReal subdiag[],PetscReal eigs[],PetscScalar V[])1144 static PetscErrorCode PetscDTSymmetricTridiagonalEigensolve(PetscInt n, PetscReal diag[], PetscReal subdiag[],
1145 PetscReal eigs[], PetscScalar V[])
1146 {
1147 char jobz = 'V'; /* eigenvalues and eigenvectors */
1148 char range = 'A'; /* all eigenvalues will be found */
1149 PetscReal VL = 0.; /* ignored because range is 'A' */
1150 PetscReal VU = 0.; /* ignored because range is 'A' */
1151 PetscBLASInt IL = 0; /* ignored because range is 'A' */
1152 PetscBLASInt IU = 0; /* ignored because range is 'A' */
1153 PetscReal abstol = 0.; /* unused */
1154 PetscBLASInt bn, bm, ldz; /* bm will equal bn on exit */
1155 PetscBLASInt *isuppz;
1156 PetscBLASInt lwork, liwork;
1157 PetscReal workquery;
1158 PetscBLASInt iworkquery;
1159 PetscBLASInt *iwork;
1160 PetscBLASInt info;
1161 PetscReal *work = NULL;
1162 PetscErrorCode ierr;
1163
1164 PetscFunctionBegin;
1165 #if !defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1166 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1167 #endif
1168 ierr = PetscBLASIntCast(n, &bn);CHKERRQ(ierr);
1169 ierr = PetscBLASIntCast(n, &ldz);CHKERRQ(ierr);
1170 #if !defined(PETSC_MISSING_LAPACK_STEGR)
1171 ierr = PetscMalloc1(2 * n, &isuppz);CHKERRQ(ierr);
1172 lwork = -1;
1173 liwork = -1;
1174 PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,&workquery,&lwork,&iworkquery,&liwork,&info));
1175 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEGR error");
1176 lwork = (PetscBLASInt) workquery;
1177 liwork = (PetscBLASInt) iworkquery;
1178 ierr = PetscMalloc2(lwork, &work, liwork, &iwork);CHKERRQ(ierr);
1179 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
1180 PetscStackCallBLAS("LAPACKstegr",LAPACKstegr_(&jobz,&range,&bn,diag,subdiag,&VL,&VU,&IL,&IU,&abstol,&bm,eigs,V,&ldz,isuppz,work,&lwork,iwork,&liwork,&info));
1181 ierr = PetscFPTrapPop();CHKERRQ(ierr);
1182 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEGR error");
1183 ierr = PetscFree2(work, iwork);CHKERRQ(ierr);
1184 ierr = PetscFree(isuppz);CHKERRQ(ierr);
1185 #elif !defined(PETSC_MISSING_LAPACK_STEQR)
1186 jobz = 'I'; /* Compute eigenvalues and eigenvectors of the
1187 tridiagonal matrix. Z is initialized to the identity
1188 matrix. */
1189 ierr = PetscMalloc1(PetscMax(1,2*n-2),&work);CHKERRQ(ierr);
1190 PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&bn,diag,subdiag,V,&ldz,work,&info));
1191 ierr = PetscFPTrapPop();CHKERRQ(ierr);
1192 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");
1193 ierr = PetscFree(work);CHKERRQ(ierr);
1194 ierr = PetscArraycpy(eigs,diag,n);CHKERRQ(ierr);
1195 #endif
1196 PetscFunctionReturn(0);
1197 }
1198
1199 /* Formula for the weights at the endpoints (-1 and 1) of Gauss-Lobatto-Jacobi
1200 * quadrature rules on the interval [-1, 1] */
PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n,PetscReal alpha,PetscReal beta,PetscReal * leftw,PetscReal * rightw)1201 static PetscErrorCode PetscDTGaussLobattoJacobiEndweights_Internal(PetscInt n, PetscReal alpha, PetscReal beta, PetscReal *leftw, PetscReal *rightw)
1202 {
1203 PetscReal twoab1;
1204 PetscInt m = n - 2;
1205 PetscReal a = alpha + 1.;
1206 PetscReal b = beta + 1.;
1207 PetscReal gra, grb;
1208
1209 PetscFunctionBegin;
1210 twoab1 = PetscPowReal(2., a + b - 1.);
1211 #if defined(PETSC_HAVE_LGAMMA)
1212 grb = PetscExpReal(2. * PetscLGamma(b+1.) + PetscLGamma(m+1.) + PetscLGamma(m+a+1.) -
1213 (PetscLGamma(m+b+1) + PetscLGamma(m+a+b+1.)));
1214 gra = PetscExpReal(2. * PetscLGamma(a+1.) + PetscLGamma(m+1.) + PetscLGamma(m+b+1.) -
1215 (PetscLGamma(m+a+1) + PetscLGamma(m+a+b+1.)));
1216 #else
1217 {
1218 PetscInt alphai = (PetscInt) alpha;
1219 PetscInt betai = (PetscInt) beta;
1220 PetscErrorCode ierr;
1221
1222 if ((PetscReal) alphai == alpha && (PetscReal) betai == beta) {
1223 PetscReal binom1, binom2;
1224
1225 ierr = PetscDTBinomial(m+b, b, &binom1);CHKERRQ(ierr);
1226 ierr = PetscDTBinomial(m+a+b, b, &binom2);CHKERRQ(ierr);
1227 grb = 1./ (binom1 * binom2);
1228 ierr = PetscDTBinomial(m+a, a, &binom1);CHKERRQ(ierr);
1229 ierr = PetscDTBinomial(m+a+b, a, &binom2);CHKERRQ(ierr);
1230 gra = 1./ (binom1 * binom2);
1231 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
1232 }
1233 #endif
1234 *leftw = twoab1 * grb / b;
1235 *rightw = twoab1 * gra / a;
1236 PetscFunctionReturn(0);
1237 }
1238
1239 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
1240 Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
PetscDTComputeJacobi(PetscReal a,PetscReal b,PetscInt n,PetscReal x,PetscReal * P)1241 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
1242 {
1243 PetscReal pn1, pn2;
1244 PetscReal cnm1, cnm1x, cnm2;
1245 PetscInt k;
1246
1247 PetscFunctionBegin;
1248 if (!n) {*P = 1.0; PetscFunctionReturn(0);}
1249 PetscDTJacobiRecurrence_Internal(1,a,b,cnm1,cnm1x,cnm2);
1250 pn2 = 1.;
1251 pn1 = cnm1 + cnm1x*x;
1252 if (n == 1) {*P = pn1; PetscFunctionReturn(0);}
1253 *P = 0.0;
1254 for (k = 2; k < n+1; ++k) {
1255 PetscDTJacobiRecurrence_Internal(k,a,b,cnm1,cnm1x,cnm2);
1256
1257 *P = (cnm1 + cnm1x*x)*pn1 - cnm2*pn2;
1258 pn2 = pn1;
1259 pn1 = *P;
1260 }
1261 PetscFunctionReturn(0);
1262 }
1263
1264 /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
PetscDTComputeJacobiDerivative(PetscReal a,PetscReal b,PetscInt n,PetscReal x,PetscInt k,PetscReal * P)1265 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscInt k, PetscReal *P)
1266 {
1267 PetscReal nP;
1268 PetscInt i;
1269 PetscErrorCode ierr;
1270
1271 PetscFunctionBegin;
1272 *P = 0.0;
1273 if (k > n) PetscFunctionReturn(0);
1274 ierr = PetscDTComputeJacobi(a+k, b+k, n-k, x, &nP);CHKERRQ(ierr);
1275 for (i = 0; i < k; i++) nP *= (a + b + n + 1. + i) * 0.5;
1276 *P = nP;
1277 PetscFunctionReturn(0);
1278 }
1279
PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints,PetscReal a,PetscReal b,PetscReal x[],PetscReal w[])1280 static PetscErrorCode PetscDTGaussJacobiQuadrature_Newton_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal x[], PetscReal w[])
1281 {
1282 PetscInt maxIter = 100;
1283 PetscReal eps = PetscExpReal(0.75 * PetscLogReal(PETSC_MACHINE_EPSILON));
1284 PetscReal a1, a6, gf;
1285 PetscInt k;
1286 PetscErrorCode ierr;
1287
1288 PetscFunctionBegin;
1289
1290 a1 = PetscPowReal(2.0, a+b+1);
1291 #if defined(PETSC_HAVE_LGAMMA)
1292 {
1293 PetscReal a2, a3, a4, a5;
1294 a2 = PetscLGamma(a + npoints + 1);
1295 a3 = PetscLGamma(b + npoints + 1);
1296 a4 = PetscLGamma(a + b + npoints + 1);
1297 a5 = PetscLGamma(npoints + 1);
1298 gf = PetscExpReal(a2 + a3 - (a4 + a5));
1299 }
1300 #else
1301 {
1302 PetscInt ia, ib;
1303
1304 ia = (PetscInt) a;
1305 ib = (PetscInt) b;
1306 gf = 1.;
1307 if (ia == a && ia >= 0) { /* compute ratio of rising factorals wrt a */
1308 for (k = 0; k < ia; k++) gf *= (npoints + 1. + k) / (npoints + b + 1. + k);
1309 } else if (b == b && ib >= 0) { /* compute ratio of rising factorials wrt b */
1310 for (k = 0; k < ib; k++) gf *= (npoints + 1. + k) / (npoints + a + 1. + k);
1311 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
1312 }
1313 #endif
1314
1315 a6 = a1 * gf;
1316 /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
1317 Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
1318 for (k = 0; k < npoints; ++k) {
1319 PetscReal r = PetscCosReal(PETSC_PI * (1. - (4.*k + 3. + 2.*b) / (4.*npoints + 2.*(a + b + 1.)))), dP;
1320 PetscInt j;
1321
1322 if (k > 0) r = 0.5 * (r + x[k-1]);
1323 for (j = 0; j < maxIter; ++j) {
1324 PetscReal s = 0.0, delta, f, fp;
1325 PetscInt i;
1326
1327 for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
1328 ierr = PetscDTComputeJacobi(a, b, npoints, r, &f);CHKERRQ(ierr);
1329 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, r, 1, &fp);CHKERRQ(ierr);
1330 delta = f / (fp - f * s);
1331 r = r - delta;
1332 if (PetscAbsReal(delta) < eps) break;
1333 }
1334 x[k] = r;
1335 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, x[k], 1, &dP);CHKERRQ(ierr);
1336 w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
1337 }
1338 PetscFunctionReturn(0);
1339 }
1340
1341 /* Compute the diagonals of the Jacobi matrix used in Golub & Welsch algorithms for Gauss-Jacobi
1342 * quadrature weight calculations on [-1,1] for exponents (1. + x)^a (1.-x)^b */
PetscDTJacobiMatrix_Internal(PetscInt nPoints,PetscReal a,PetscReal b,PetscReal * d,PetscReal * s)1343 static PetscErrorCode PetscDTJacobiMatrix_Internal(PetscInt nPoints, PetscReal a, PetscReal b, PetscReal *d, PetscReal *s)
1344 {
1345 PetscInt i;
1346
1347 PetscFunctionBegin;
1348 for (i = 0; i < nPoints; i++) {
1349 PetscReal A, B, C;
1350
1351 PetscDTJacobiRecurrence_Internal(i+1,a,b,A,B,C);
1352 d[i] = -A / B;
1353 if (i) s[i-1] *= C / B;
1354 if (i < nPoints - 1) s[i] = 1. / B;
1355 }
1356 PetscFunctionReturn(0);
1357 }
1358
PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints,PetscReal a,PetscReal b,PetscReal * x,PetscReal * w)1359 static PetscErrorCode PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
1360 {
1361 PetscReal mu0;
1362 PetscReal ga, gb, gab;
1363 PetscInt i;
1364 PetscErrorCode ierr;
1365
1366 PetscFunctionBegin;
1367 ierr = PetscCitationsRegister(GolubWelschCitation, &GolubWelschCite);CHKERRQ(ierr);
1368
1369 #if defined(PETSC_HAVE_TGAMMA)
1370 ga = PetscTGamma(a + 1);
1371 gb = PetscTGamma(b + 1);
1372 gab = PetscTGamma(a + b + 2);
1373 #else
1374 {
1375 PetscInt ia, ib;
1376
1377 ia = (PetscInt) a;
1378 ib = (PetscInt) b;
1379 if (ia == a && ib == b && ia + 1 > 0 && ib + 1 > 0 && ia + ib + 2 > 0) { /* All gamma(x) terms are (x-1)! terms */
1380 ierr = PetscDTFactorial(ia, &ga);CHKERRQ(ierr);
1381 ierr = PetscDTFactorial(ib, &gb);CHKERRQ(ierr);
1382 ierr = PetscDTFactorial(ia + ib + 1, &gb);CHKERRQ(ierr);
1383 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
1384 }
1385 #endif
1386 mu0 = PetscPowReal(2.,a + b + 1.) * ga * gb / gab;
1387
1388 #if defined(PETSCDTGAUSSIANQUADRATURE_EIG)
1389 {
1390 PetscReal *diag, *subdiag;
1391 PetscScalar *V;
1392
1393 ierr = PetscMalloc2(npoints, &diag, npoints, &subdiag);CHKERRQ(ierr);
1394 ierr = PetscMalloc1(npoints*npoints, &V);CHKERRQ(ierr);
1395 ierr = PetscDTJacobiMatrix_Internal(npoints, a, b, diag, subdiag);CHKERRQ(ierr);
1396 for (i = 0; i < npoints - 1; i++) subdiag[i] = PetscSqrtReal(subdiag[i]);
1397 ierr = PetscDTSymmetricTridiagonalEigensolve(npoints, diag, subdiag, x, V);CHKERRQ(ierr);
1398 for (i = 0; i < npoints; i++) w[i] = PetscSqr(PetscRealPart(V[i * npoints])) * mu0;
1399 ierr = PetscFree(V);CHKERRQ(ierr);
1400 ierr = PetscFree2(diag, subdiag);CHKERRQ(ierr);
1401 }
1402 #else
1403 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "A LAPACK symmetric tridiagonal eigensolver could not be found");
1404 #endif
1405 { /* As of March 2, 2020, The Sun Performance Library breaks the LAPACK contract for xstegr and xsteqr: the
1406 eigenvalues are not guaranteed to be in ascending order. So we heave a passive aggressive sigh and check that
1407 the eigenvalues are sorted */
1408 PetscBool sorted;
1409
1410 ierr = PetscSortedReal(npoints, x, &sorted);CHKERRQ(ierr);
1411 if (!sorted) {
1412 PetscInt *order, i;
1413 PetscReal *tmp;
1414
1415 ierr = PetscMalloc2(npoints, &order, npoints, &tmp);CHKERRQ(ierr);
1416 for (i = 0; i < npoints; i++) order[i] = i;
1417 ierr = PetscSortRealWithPermutation(npoints, x, order);CHKERRQ(ierr);
1418 ierr = PetscArraycpy(tmp, x, npoints);CHKERRQ(ierr);
1419 for (i = 0; i < npoints; i++) x[i] = tmp[order[i]];
1420 ierr = PetscArraycpy(tmp, w, npoints);CHKERRQ(ierr);
1421 for (i = 0; i < npoints; i++) w[i] = tmp[order[i]];
1422 ierr = PetscFree2(order, tmp);CHKERRQ(ierr);
1423 }
1424 }
1425 PetscFunctionReturn(0);
1426 }
1427
PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha,PetscReal beta,PetscReal x[],PetscReal w[],PetscBool newton)1428 static PetscErrorCode PetscDTGaussJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1429 {
1430 PetscErrorCode ierr;
1431
1432 PetscFunctionBegin;
1433 if (npoints < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be positive");
1434 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1435 if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1.");
1436 if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1.");
1437
1438 if (newton) {
1439 ierr = PetscDTGaussJacobiQuadrature_Newton_Internal(npoints, alpha, beta, x, w);CHKERRQ(ierr);
1440 } else {
1441 ierr = PetscDTGaussJacobiQuadrature_GolubWelsch_Internal(npoints, alpha, beta, x, w);CHKERRQ(ierr);
1442 }
1443 if (alpha == beta) { /* symmetrize */
1444 PetscInt i;
1445 for (i = 0; i < (npoints + 1) / 2; i++) {
1446 PetscInt j = npoints - 1 - i;
1447 PetscReal xi = x[i];
1448 PetscReal xj = x[j];
1449 PetscReal wi = w[i];
1450 PetscReal wj = w[j];
1451
1452 x[i] = (xi - xj) / 2.;
1453 x[j] = (xj - xi) / 2.;
1454 w[i] = w[j] = (wi + wj) / 2.;
1455 }
1456 }
1457 PetscFunctionReturn(0);
1458 }
1459
1460 /*@
1461 PetscDTGaussJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1462 $(x-a)^\alpha (x-b)^\beta$.
1463
1464 Not collective
1465
1466 Input Parameters:
1467 + npoints - the number of points in the quadrature rule
1468 . a - the left endpoint of the interval
1469 . b - the right endpoint of the interval
1470 . alpha - the left exponent
1471 - beta - the right exponent
1472
1473 Output Parameters:
1474 + x - array of length npoints, the locations of the quadrature points
1475 - w - array of length npoints, the weights of the quadrature points
1476
1477 Level: intermediate
1478
1479 Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 1.
1480 @*/
PetscDTGaussJacobiQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal alpha,PetscReal beta,PetscReal x[],PetscReal w[])1481 PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1482 {
1483 PetscInt i;
1484 PetscErrorCode ierr;
1485
1486 PetscFunctionBegin;
1487 ierr = PetscDTGaussJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);CHKERRQ(ierr);
1488 if (a != -1. || b != 1.) { /* shift */
1489 for (i = 0; i < npoints; i++) {
1490 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1491 w[i] *= (b - a) / 2.;
1492 }
1493 }
1494 PetscFunctionReturn(0);
1495 }
1496
PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha,PetscReal beta,PetscReal x[],PetscReal w[],PetscBool newton)1497 static PetscErrorCode PetscDTGaussLobattoJacobiQuadrature_Internal(PetscInt npoints,PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[], PetscBool newton)
1498 {
1499 PetscInt i;
1500 PetscErrorCode ierr;
1501
1502 PetscFunctionBegin;
1503 if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be positive");
1504 /* If asking for a 1D Lobatto point, just return the non-Lobatto 1D point */
1505 if (alpha <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"alpha must be > -1.");
1506 if (beta <= -1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"beta must be > -1.");
1507
1508 x[0] = -1.;
1509 x[npoints-1] = 1.;
1510 if (npoints > 2) {
1511 ierr = PetscDTGaussJacobiQuadrature_Internal(npoints-2, alpha+1., beta+1., &x[1], &w[1], newton);CHKERRQ(ierr);
1512 }
1513 for (i = 1; i < npoints - 1; i++) {
1514 w[i] /= (1. - x[i]*x[i]);
1515 }
1516 ierr = PetscDTGaussLobattoJacobiEndweights_Internal(npoints, alpha, beta, &w[0], &w[npoints-1]);CHKERRQ(ierr);
1517 PetscFunctionReturn(0);
1518 }
1519
1520 /*@
1521 PetscDTGaussLobattoJacobiQuadrature - quadrature for the interval [a, b] with the weight function
1522 $(x-a)^\alpha (x-b)^\beta$, with endpoints a and b included as quadrature points.
1523
1524 Not collective
1525
1526 Input Parameters:
1527 + npoints - the number of points in the quadrature rule
1528 . a - the left endpoint of the interval
1529 . b - the right endpoint of the interval
1530 . alpha - the left exponent
1531 - beta - the right exponent
1532
1533 Output Parameters:
1534 + x - array of length npoints, the locations of the quadrature points
1535 - w - array of length npoints, the weights of the quadrature points
1536
1537 Level: intermediate
1538
1539 Note: this quadrature rule is exact for polynomials up to degree 2*npoints - 3.
1540 @*/
PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal alpha,PetscReal beta,PetscReal x[],PetscReal w[])1541 PetscErrorCode PetscDTGaussLobattoJacobiQuadrature(PetscInt npoints,PetscReal a, PetscReal b, PetscReal alpha, PetscReal beta, PetscReal x[], PetscReal w[])
1542 {
1543 PetscInt i;
1544 PetscErrorCode ierr;
1545
1546 PetscFunctionBegin;
1547 ierr = PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, alpha, beta, x, w, PetscDTGaussQuadratureNewton_Internal);CHKERRQ(ierr);
1548 if (a != -1. || b != 1.) { /* shift */
1549 for (i = 0; i < npoints; i++) {
1550 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1551 w[i] *= (b - a) / 2.;
1552 }
1553 }
1554 PetscFunctionReturn(0);
1555 }
1556
1557 /*@
1558 PetscDTGaussQuadrature - create Gauss-Legendre quadrature
1559
1560 Not Collective
1561
1562 Input Arguments:
1563 + npoints - number of points
1564 . a - left end of interval (often-1)
1565 - b - right end of interval (often +1)
1566
1567 Output Arguments:
1568 + x - quadrature points
1569 - w - quadrature weights
1570
1571 Level: intermediate
1572
1573 References:
1574 . 1. - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.
1575
1576 .seealso: PetscDTLegendreEval()
1577 @*/
PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal * x,PetscReal * w)1578 PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
1579 {
1580 PetscInt i;
1581 PetscErrorCode ierr;
1582
1583 PetscFunctionBegin;
1584 ierr = PetscDTGaussJacobiQuadrature_Internal(npoints, 0., 0., x, w, PetscDTGaussQuadratureNewton_Internal);CHKERRQ(ierr);
1585 if (a != -1. || b != 1.) { /* shift */
1586 for (i = 0; i < npoints; i++) {
1587 x[i] = (x[i] + 1.) * ((b - a) / 2.) + a;
1588 w[i] *= (b - a) / 2.;
1589 }
1590 }
1591 PetscFunctionReturn(0);
1592 }
1593
1594 /*@C
1595 PetscDTGaussLobattoLegendreQuadrature - creates a set of the locations and weights of the Gauss-Lobatto-Legendre
1596 nodes of a given size on the domain [-1,1]
1597
1598 Not Collective
1599
1600 Input Parameter:
1601 + n - number of grid nodes
1602 - type - PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA or PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON
1603
1604 Output Arguments:
1605 + x - quadrature points
1606 - w - quadrature weights
1607
1608 Notes:
1609 For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
1610 close enough to the desired solution
1611
1612 These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes
1613
1614 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
1615
1616 Level: intermediate
1617
1618 .seealso: PetscDTGaussQuadrature()
1619
1620 @*/
PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal * x,PetscReal * w)1621 PetscErrorCode PetscDTGaussLobattoLegendreQuadrature(PetscInt npoints,PetscGaussLobattoLegendreCreateType type,PetscReal *x,PetscReal *w)
1622 {
1623 PetscBool newton;
1624 PetscErrorCode ierr;
1625
1626 PetscFunctionBegin;
1627 if (npoints < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide at least 2 grid points per element");
1628 newton = (PetscBool) (type == PETSCGAUSSLOBATTOLEGENDRE_VIA_NEWTON);
1629 ierr = PetscDTGaussLobattoJacobiQuadrature_Internal(npoints, 0., 0., x, w, newton);CHKERRQ(ierr);
1630 PetscFunctionReturn(0);
1631 }
1632
1633 /*@
1634 PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
1635
1636 Not Collective
1637
1638 Input Arguments:
1639 + dim - The spatial dimension
1640 . Nc - The number of components
1641 . npoints - number of points in one dimension
1642 . a - left end of interval (often-1)
1643 - b - right end of interval (often +1)
1644
1645 Output Argument:
1646 . q - A PetscQuadrature object
1647
1648 Level: intermediate
1649
1650 .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
1651 @*/
PetscDTGaussTensorQuadrature(PetscInt dim,PetscInt Nc,PetscInt npoints,PetscReal a,PetscReal b,PetscQuadrature * q)1652 PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1653 {
1654 PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
1655 PetscReal *x, *w, *xw, *ww;
1656 PetscErrorCode ierr;
1657
1658 PetscFunctionBegin;
1659 ierr = PetscMalloc1(totpoints*dim,&x);CHKERRQ(ierr);
1660 ierr = PetscMalloc1(totpoints*Nc,&w);CHKERRQ(ierr);
1661 /* Set up the Golub-Welsch system */
1662 switch (dim) {
1663 case 0:
1664 ierr = PetscFree(x);CHKERRQ(ierr);
1665 ierr = PetscFree(w);CHKERRQ(ierr);
1666 ierr = PetscMalloc1(1, &x);CHKERRQ(ierr);
1667 ierr = PetscMalloc1(Nc, &w);CHKERRQ(ierr);
1668 x[0] = 0.0;
1669 for (c = 0; c < Nc; ++c) w[c] = 1.0;
1670 break;
1671 case 1:
1672 ierr = PetscMalloc1(npoints,&ww);CHKERRQ(ierr);
1673 ierr = PetscDTGaussQuadrature(npoints, a, b, x, ww);CHKERRQ(ierr);
1674 for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
1675 ierr = PetscFree(ww);CHKERRQ(ierr);
1676 break;
1677 case 2:
1678 ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr);
1679 ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr);
1680 for (i = 0; i < npoints; ++i) {
1681 for (j = 0; j < npoints; ++j) {
1682 x[(i*npoints+j)*dim+0] = xw[i];
1683 x[(i*npoints+j)*dim+1] = xw[j];
1684 for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
1685 }
1686 }
1687 ierr = PetscFree2(xw,ww);CHKERRQ(ierr);
1688 break;
1689 case 3:
1690 ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr);
1691 ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr);
1692 for (i = 0; i < npoints; ++i) {
1693 for (j = 0; j < npoints; ++j) {
1694 for (k = 0; k < npoints; ++k) {
1695 x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
1696 x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
1697 x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
1698 for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
1699 }
1700 }
1701 }
1702 ierr = PetscFree2(xw,ww);CHKERRQ(ierr);
1703 break;
1704 default:
1705 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
1706 }
1707 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr);
1708 ierr = PetscQuadratureSetOrder(*q, 2*npoints-1);CHKERRQ(ierr);
1709 ierr = PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);CHKERRQ(ierr);
1710 ierr = PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");CHKERRQ(ierr);
1711 PetscFunctionReturn(0);
1712 }
1713
1714 /*@
1715 PetscDTStroudConicalQuadrature - create Stroud conical quadrature for a simplex
1716
1717 Not Collective
1718
1719 Input Arguments:
1720 + dim - The simplex dimension
1721 . Nc - The number of components
1722 . npoints - The number of points in one dimension
1723 . a - left end of interval (often-1)
1724 - b - right end of interval (often +1)
1725
1726 Output Argument:
1727 . q - A PetscQuadrature object
1728
1729 Level: intermediate
1730
1731 References:
1732 . 1. - Karniadakis and Sherwin. FIAT
1733
1734 Note: For dim == 1, this is Gauss-Legendre quadrature
1735
1736 .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
1737 @*/
PetscDTStroudConicalQuadrature(PetscInt dim,PetscInt Nc,PetscInt npoints,PetscReal a,PetscReal b,PetscQuadrature * q)1738 PetscErrorCode PetscDTStroudConicalQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
1739 {
1740 PetscInt totprev, totrem;
1741 PetscInt totpoints;
1742 PetscReal *p1, *w1;
1743 PetscReal *x, *w;
1744 PetscInt i, j, k, l, m, pt, c;
1745 PetscErrorCode ierr;
1746
1747 PetscFunctionBegin;
1748 if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
1749 totpoints = 1;
1750 for (i = 0, totpoints = 1; i < dim; i++) totpoints *= npoints;
1751 ierr = PetscMalloc1(totpoints*dim, &x);CHKERRQ(ierr);
1752 ierr = PetscMalloc1(totpoints*Nc, &w);CHKERRQ(ierr);
1753 ierr = PetscMalloc2(npoints, &p1, npoints, &w1);CHKERRQ(ierr);
1754 for (i = 0; i < totpoints*Nc; i++) w[i] = 1.;
1755 for (i = 0, totprev = 1, totrem = totpoints / npoints; i < dim; i++) {
1756 PetscReal mul;
1757
1758 mul = PetscPowReal(2.,-i);
1759 ierr = PetscDTGaussJacobiQuadrature(npoints, -1., 1., i, 0.0, p1, w1);CHKERRQ(ierr);
1760 for (pt = 0, l = 0; l < totprev; l++) {
1761 for (j = 0; j < npoints; j++) {
1762 for (m = 0; m < totrem; m++, pt++) {
1763 for (k = 0; k < i; k++) x[pt*dim+k] = (x[pt*dim+k]+1.)*(1.-p1[j])*0.5 - 1.;
1764 x[pt * dim + i] = p1[j];
1765 for (c = 0; c < Nc; c++) w[pt*Nc + c] *= mul * w1[j];
1766 }
1767 }
1768 }
1769 totprev *= npoints;
1770 totrem /= npoints;
1771 }
1772 ierr = PetscFree2(p1, w1);CHKERRQ(ierr);
1773 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr);
1774 ierr = PetscQuadratureSetOrder(*q, 2*npoints-1);CHKERRQ(ierr);
1775 ierr = PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);CHKERRQ(ierr);
1776 ierr = PetscObjectChangeTypeName((PetscObject)*q,"StroudConical");CHKERRQ(ierr);
1777 PetscFunctionReturn(0);
1778 }
1779
1780 /*@
1781 PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell
1782
1783 Not Collective
1784
1785 Input Arguments:
1786 + dim - The cell dimension
1787 . level - The number of points in one dimension, 2^l
1788 . a - left end of interval (often-1)
1789 - b - right end of interval (often +1)
1790
1791 Output Argument:
1792 . q - A PetscQuadrature object
1793
1794 Level: intermediate
1795
1796 .seealso: PetscDTGaussTensorQuadrature()
1797 @*/
PetscDTTanhSinhTensorQuadrature(PetscInt dim,PetscInt level,PetscReal a,PetscReal b,PetscQuadrature * q)1798 PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
1799 {
1800 const PetscInt p = 16; /* Digits of precision in the evaluation */
1801 const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
1802 const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
1803 const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */
1804 PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */
1805 PetscReal wk = 0.5*PETSC_PI; /* Quadrature weight at x_k */
1806 PetscReal *x, *w;
1807 PetscInt K, k, npoints;
1808 PetscErrorCode ierr;
1809
1810 PetscFunctionBegin;
1811 if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
1812 if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
1813 /* Find K such that the weights are < 32 digits of precision */
1814 for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
1815 wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
1816 }
1817 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr);
1818 ierr = PetscQuadratureSetOrder(*q, 2*K+1);CHKERRQ(ierr);
1819 npoints = 2*K-1;
1820 ierr = PetscMalloc1(npoints*dim, &x);CHKERRQ(ierr);
1821 ierr = PetscMalloc1(npoints, &w);CHKERRQ(ierr);
1822 /* Center term */
1823 x[0] = beta;
1824 w[0] = 0.5*alpha*PETSC_PI;
1825 for (k = 1; k < K; ++k) {
1826 wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1827 xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
1828 x[2*k-1] = -alpha*xk+beta;
1829 w[2*k-1] = wk;
1830 x[2*k+0] = alpha*xk+beta;
1831 w[2*k+0] = wk;
1832 }
1833 ierr = PetscQuadratureSetData(*q, dim, 1, npoints, x, w);CHKERRQ(ierr);
1834 PetscFunctionReturn(0);
1835 }
1836
PetscDTTanhSinhIntegrate(void (* func)(PetscReal,PetscReal *),PetscReal a,PetscReal b,PetscInt digits,PetscReal * sol)1837 PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1838 {
1839 const PetscInt p = 16; /* Digits of precision in the evaluation */
1840 const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
1841 const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
1842 PetscReal h = 1.0; /* Step size, length between x_k */
1843 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
1844 PetscReal osum = 0.0; /* Integral on last level */
1845 PetscReal psum = 0.0; /* Integral on the level before the last level */
1846 PetscReal sum; /* Integral on current level */
1847 PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1848 PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1849 PetscReal wk; /* Quadrature weight at x_k */
1850 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
1851 PetscInt d; /* Digits of precision in the integral */
1852
1853 PetscFunctionBegin;
1854 if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1855 /* Center term */
1856 func(beta, &lval);
1857 sum = 0.5*alpha*PETSC_PI*lval;
1858 /* */
1859 do {
1860 PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
1861 PetscInt k = 1;
1862
1863 ++l;
1864 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1865 /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1866 psum = osum;
1867 osum = sum;
1868 h *= 0.5;
1869 sum *= 0.5;
1870 do {
1871 wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1872 yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
1873 lx = -alpha*(1.0 - yk)+beta;
1874 rx = alpha*(1.0 - yk)+beta;
1875 func(lx, &lval);
1876 func(rx, &rval);
1877 lterm = alpha*wk*lval;
1878 maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
1879 sum += lterm;
1880 rterm = alpha*wk*rval;
1881 maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
1882 sum += rterm;
1883 ++k;
1884 /* Only need to evaluate every other point on refined levels */
1885 if (l != 1) ++k;
1886 } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */
1887
1888 d1 = PetscLog10Real(PetscAbsReal(sum - osum));
1889 d2 = PetscLog10Real(PetscAbsReal(sum - psum));
1890 d3 = PetscLog10Real(maxTerm) - p;
1891 if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
1892 else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
1893 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1894 } while (d < digits && l < 12);
1895 *sol = sum;
1896
1897 PetscFunctionReturn(0);
1898 }
1899
1900 #if defined(PETSC_HAVE_MPFR)
PetscDTTanhSinhIntegrateMPFR(void (* func)(PetscReal,PetscReal *),PetscReal a,PetscReal b,PetscInt digits,PetscReal * sol)1901 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1902 {
1903 const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */
1904 PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
1905 mpfr_t alpha; /* Half-width of the integration interval */
1906 mpfr_t beta; /* Center of the integration interval */
1907 mpfr_t h; /* Step size, length between x_k */
1908 mpfr_t osum; /* Integral on last level */
1909 mpfr_t psum; /* Integral on the level before the last level */
1910 mpfr_t sum; /* Integral on current level */
1911 mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
1912 mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
1913 mpfr_t wk; /* Quadrature weight at x_k */
1914 PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
1915 PetscInt d; /* Digits of precision in the integral */
1916 mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;
1917
1918 PetscFunctionBegin;
1919 if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
1920 /* Create high precision storage */
1921 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);
1922 /* Initialization */
1923 mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
1924 mpfr_set_d(beta, 0.5*(b+a), MPFR_RNDN);
1925 mpfr_set_d(osum, 0.0, MPFR_RNDN);
1926 mpfr_set_d(psum, 0.0, MPFR_RNDN);
1927 mpfr_set_d(h, 1.0, MPFR_RNDN);
1928 mpfr_const_pi(pi2, MPFR_RNDN);
1929 mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
1930 /* Center term */
1931 func(0.5*(b+a), &lval);
1932 mpfr_set(sum, pi2, MPFR_RNDN);
1933 mpfr_mul(sum, sum, alpha, MPFR_RNDN);
1934 mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
1935 /* */
1936 do {
1937 PetscReal d1, d2, d3, d4;
1938 PetscInt k = 1;
1939
1940 ++l;
1941 mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
1942 /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1943 /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1944 mpfr_set(psum, osum, MPFR_RNDN);
1945 mpfr_set(osum, sum, MPFR_RNDN);
1946 mpfr_mul_d(h, h, 0.5, MPFR_RNDN);
1947 mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
1948 do {
1949 mpfr_set_si(kh, k, MPFR_RNDN);
1950 mpfr_mul(kh, kh, h, MPFR_RNDN);
1951 /* Weight */
1952 mpfr_set(wk, h, MPFR_RNDN);
1953 mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
1954 mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
1955 mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
1956 mpfr_cosh(tmp, msinh, MPFR_RNDN);
1957 mpfr_sqr(tmp, tmp, MPFR_RNDN);
1958 mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
1959 mpfr_div(wk, wk, tmp, MPFR_RNDN);
1960 /* Abscissa */
1961 mpfr_set_d(yk, 1.0, MPFR_RNDZ);
1962 mpfr_cosh(tmp, msinh, MPFR_RNDN);
1963 mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1964 mpfr_exp(tmp, msinh, MPFR_RNDN);
1965 mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1966 /* Quadrature points */
1967 mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
1968 mpfr_mul(lx, lx, alpha, MPFR_RNDU);
1969 mpfr_add(lx, lx, beta, MPFR_RNDU);
1970 mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
1971 mpfr_mul(rx, rx, alpha, MPFR_RNDD);
1972 mpfr_add(rx, rx, beta, MPFR_RNDD);
1973 /* Evaluation */
1974 func(mpfr_get_d(lx, MPFR_RNDU), &lval);
1975 func(mpfr_get_d(rx, MPFR_RNDD), &rval);
1976 /* Update */
1977 mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1978 mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
1979 mpfr_add(sum, sum, tmp, MPFR_RNDN);
1980 mpfr_abs(tmp, tmp, MPFR_RNDN);
1981 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1982 mpfr_set(curTerm, tmp, MPFR_RNDN);
1983 mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1984 mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
1985 mpfr_add(sum, sum, tmp, MPFR_RNDN);
1986 mpfr_abs(tmp, tmp, MPFR_RNDN);
1987 mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1988 mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
1989 ++k;
1990 /* Only need to evaluate every other point on refined levels */
1991 if (l != 1) ++k;
1992 mpfr_log10(tmp, wk, MPFR_RNDN);
1993 mpfr_abs(tmp, tmp, MPFR_RNDN);
1994 } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
1995 mpfr_sub(tmp, sum, osum, MPFR_RNDN);
1996 mpfr_abs(tmp, tmp, MPFR_RNDN);
1997 mpfr_log10(tmp, tmp, MPFR_RNDN);
1998 d1 = mpfr_get_d(tmp, MPFR_RNDN);
1999 mpfr_sub(tmp, sum, psum, MPFR_RNDN);
2000 mpfr_abs(tmp, tmp, MPFR_RNDN);
2001 mpfr_log10(tmp, tmp, MPFR_RNDN);
2002 d2 = mpfr_get_d(tmp, MPFR_RNDN);
2003 mpfr_log10(tmp, maxTerm, MPFR_RNDN);
2004 d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
2005 mpfr_log10(tmp, curTerm, MPFR_RNDN);
2006 d4 = mpfr_get_d(tmp, MPFR_RNDN);
2007 d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
2008 } while (d < digits && l < 8);
2009 *sol = mpfr_get_d(sum, MPFR_RNDN);
2010 /* Cleanup */
2011 mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
2012 PetscFunctionReturn(0);
2013 }
2014 #else
2015
PetscDTTanhSinhIntegrateMPFR(void (* func)(PetscReal,PetscReal *),PetscReal a,PetscReal b,PetscInt digits,PetscReal * sol)2016 PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
2017 {
2018 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
2019 }
2020 #endif
2021
2022 /* Overwrites A. Can only handle full-rank problems with m>=n
2023 * A in column-major format
2024 * Ainv in row-major format
2025 * tau has length m
2026 * worksize must be >= max(1,n)
2027 */
PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal * A_in,PetscReal * Ainv_out,PetscScalar * tau,PetscInt worksize,PetscScalar * work)2028 static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2029 {
2030 PetscErrorCode ierr;
2031 PetscBLASInt M,N,K,lda,ldb,ldwork,info;
2032 PetscScalar *A,*Ainv,*R,*Q,Alpha;
2033
2034 PetscFunctionBegin;
2035 #if defined(PETSC_USE_COMPLEX)
2036 {
2037 PetscInt i,j;
2038 ierr = PetscMalloc2(m*n,&A,m*n,&Ainv);CHKERRQ(ierr);
2039 for (j=0; j<n; j++) {
2040 for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
2041 }
2042 mstride = m;
2043 }
2044 #else
2045 A = A_in;
2046 Ainv = Ainv_out;
2047 #endif
2048
2049 ierr = PetscBLASIntCast(m,&M);CHKERRQ(ierr);
2050 ierr = PetscBLASIntCast(n,&N);CHKERRQ(ierr);
2051 ierr = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr);
2052 ierr = PetscBLASIntCast(worksize,&ldwork);CHKERRQ(ierr);
2053 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
2054 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
2055 ierr = PetscFPTrapPop();CHKERRQ(ierr);
2056 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2057 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2058
2059 /* Extract an explicit representation of Q */
2060 Q = Ainv;
2061 ierr = PetscArraycpy(Q,A,mstride*n);CHKERRQ(ierr);
2062 K = N; /* full rank */
2063 PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
2064 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
2065
2066 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2067 Alpha = 1.0;
2068 ldb = lda;
2069 PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
2070 /* Ainv is Q, overwritten with inverse */
2071
2072 #if defined(PETSC_USE_COMPLEX)
2073 {
2074 PetscInt i;
2075 for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
2076 ierr = PetscFree2(A,Ainv);CHKERRQ(ierr);
2077 }
2078 #endif
2079 PetscFunctionReturn(0);
2080 }
2081
2082 /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal * x,PetscInt ndegree,const PetscInt * degrees,PetscBool Transpose,PetscReal * B)2083 static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
2084 {
2085 PetscErrorCode ierr;
2086 PetscReal *Bv;
2087 PetscInt i,j;
2088
2089 PetscFunctionBegin;
2090 ierr = PetscMalloc1((ninterval+1)*ndegree,&Bv);CHKERRQ(ierr);
2091 /* Point evaluation of L_p on all the source vertices */
2092 ierr = PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);CHKERRQ(ierr);
2093 /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
2094 for (i=0; i<ninterval; i++) {
2095 for (j=0; j<ndegree; j++) {
2096 if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
2097 else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
2098 }
2099 }
2100 ierr = PetscFree(Bv);CHKERRQ(ierr);
2101 PetscFunctionReturn(0);
2102 }
2103
2104 /*@
2105 PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
2106
2107 Not Collective
2108
2109 Input Arguments:
2110 + degree - degree of reconstruction polynomial
2111 . nsource - number of source intervals
2112 . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
2113 . ntarget - number of target intervals
2114 - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
2115
2116 Output Arguments:
2117 . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
2118
2119 Level: advanced
2120
2121 .seealso: PetscDTLegendreEval()
2122 @*/
PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal * sourcex,PetscInt ntarget,const PetscReal * targetx,PetscReal * R)2123 PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
2124 {
2125 PetscErrorCode ierr;
2126 PetscInt i,j,k,*bdegrees,worksize;
2127 PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
2128 PetscScalar *tau,*work;
2129
2130 PetscFunctionBegin;
2131 PetscValidRealPointer(sourcex,3);
2132 PetscValidRealPointer(targetx,5);
2133 PetscValidRealPointer(R,6);
2134 if (degree >= nsource) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Reconstruction degree %D must be less than number of source intervals %D",degree,nsource);
2135 if (PetscDefined(USE_DEBUG)) {
2136 for (i=0; i<nsource; i++) {
2137 if (sourcex[i] >= sourcex[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Source interval %D has negative orientation (%g,%g)",i,(double)sourcex[i],(double)sourcex[i+1]);
2138 }
2139 for (i=0; i<ntarget; i++) {
2140 if (targetx[i] >= targetx[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Target interval %D has negative orientation (%g,%g)",i,(double)targetx[i],(double)targetx[i+1]);
2141 }
2142 }
2143 xmin = PetscMin(sourcex[0],targetx[0]);
2144 xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
2145 center = (xmin + xmax)/2;
2146 hscale = (xmax - xmin)/2;
2147 worksize = nsource;
2148 ierr = PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);CHKERRQ(ierr);
2149 ierr = PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);CHKERRQ(ierr);
2150 for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
2151 for (i=0; i<=degree; i++) bdegrees[i] = i+1;
2152 ierr = PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);CHKERRQ(ierr);
2153 ierr = PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);CHKERRQ(ierr);
2154 for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
2155 ierr = PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);CHKERRQ(ierr);
2156 for (i=0; i<ntarget; i++) {
2157 PetscReal rowsum = 0;
2158 for (j=0; j<nsource; j++) {
2159 PetscReal sum = 0;
2160 for (k=0; k<degree+1; k++) {
2161 sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
2162 }
2163 R[i*nsource+j] = sum;
2164 rowsum += sum;
2165 }
2166 for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
2167 }
2168 ierr = PetscFree4(bdegrees,sourcey,Bsource,work);CHKERRQ(ierr);
2169 ierr = PetscFree4(tau,Bsinv,targety,Btarget);CHKERRQ(ierr);
2170 PetscFunctionReturn(0);
2171 }
2172
2173 /*@C
2174 PetscGaussLobattoLegendreIntegrate - Compute the L2 integral of a function on the GLL points
2175
2176 Not Collective
2177
2178 Input Parameter:
2179 + n - the number of GLL nodes
2180 . nodes - the GLL nodes
2181 . weights - the GLL weights
2182 - f - the function values at the nodes
2183
2184 Output Parameter:
2185 . in - the value of the integral
2186
2187 Level: beginner
2188
2189 .seealso: PetscDTGaussLobattoLegendreQuadrature()
2190
2191 @*/
PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal * nodes,PetscReal * weights,const PetscReal * f,PetscReal * in)2192 PetscErrorCode PetscGaussLobattoLegendreIntegrate(PetscInt n,PetscReal *nodes,PetscReal *weights,const PetscReal *f,PetscReal *in)
2193 {
2194 PetscInt i;
2195
2196 PetscFunctionBegin;
2197 *in = 0.;
2198 for (i=0; i<n; i++) {
2199 *in += f[i]*f[i]*weights[i];
2200 }
2201 PetscFunctionReturn(0);
2202 }
2203
2204 /*@C
2205 PetscGaussLobattoLegendreElementLaplacianCreate - computes the Laplacian for a single 1d GLL element
2206
2207 Not Collective
2208
2209 Input Parameter:
2210 + n - the number of GLL nodes
2211 . nodes - the GLL nodes
2212 - weights - the GLL weights
2213
2214 Output Parameter:
2215 . A - the stiffness element
2216
2217 Level: beginner
2218
2219 Notes:
2220 Destroy this with PetscGaussLobattoLegendreElementLaplacianDestroy()
2221
2222 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)
2223
2224 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()
2225
2226 @*/
PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal * nodes,PetscReal * weights,PetscReal *** AA)2227 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2228 {
2229 PetscReal **A;
2230 PetscErrorCode ierr;
2231 const PetscReal *gllnodes = nodes;
2232 const PetscInt p = n-1;
2233 PetscReal z0,z1,z2 = -1,x,Lpj,Lpr;
2234 PetscInt i,j,nn,r;
2235
2236 PetscFunctionBegin;
2237 ierr = PetscMalloc1(n,&A);CHKERRQ(ierr);
2238 ierr = PetscMalloc1(n*n,&A[0]);CHKERRQ(ierr);
2239 for (i=1; i<n; i++) A[i] = A[i-1]+n;
2240
2241 for (j=1; j<p; j++) {
2242 x = gllnodes[j];
2243 z0 = 1.;
2244 z1 = x;
2245 for (nn=1; nn<p; nn++) {
2246 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2247 z0 = z1;
2248 z1 = z2;
2249 }
2250 Lpj=z2;
2251 for (r=1; r<p; r++) {
2252 if (r == j) {
2253 A[j][j]=2./(3.*(1.-gllnodes[j]*gllnodes[j])*Lpj*Lpj);
2254 } else {
2255 x = gllnodes[r];
2256 z0 = 1.;
2257 z1 = x;
2258 for (nn=1; nn<p; nn++) {
2259 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2260 z0 = z1;
2261 z1 = z2;
2262 }
2263 Lpr = z2;
2264 A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(gllnodes[j]-gllnodes[r])*(gllnodes[j]-gllnodes[r]));
2265 }
2266 }
2267 }
2268 for (j=1; j<p+1; j++) {
2269 x = gllnodes[j];
2270 z0 = 1.;
2271 z1 = x;
2272 for (nn=1; nn<p; nn++) {
2273 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2274 z0 = z1;
2275 z1 = z2;
2276 }
2277 Lpj = z2;
2278 A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+gllnodes[j])*(1.+gllnodes[j]));
2279 A[0][j] = A[j][0];
2280 }
2281 for (j=0; j<p; j++) {
2282 x = gllnodes[j];
2283 z0 = 1.;
2284 z1 = x;
2285 for (nn=1; nn<p; nn++) {
2286 z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
2287 z0 = z1;
2288 z1 = z2;
2289 }
2290 Lpj=z2;
2291
2292 A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-gllnodes[j])*(1.-gllnodes[j]));
2293 A[j][p] = A[p][j];
2294 }
2295 A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.;
2296 A[p][p]=A[0][0];
2297 *AA = A;
2298 PetscFunctionReturn(0);
2299 }
2300
2301 /*@C
2302 PetscGaussLobattoLegendreElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element
2303
2304 Not Collective
2305
2306 Input Parameter:
2307 + n - the number of GLL nodes
2308 . nodes - the GLL nodes
2309 . weights - the GLL weightss
2310 - A - the stiffness element
2311
2312 Level: beginner
2313
2314 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate()
2315
2316 @*/
PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal * nodes,PetscReal * weights,PetscReal *** AA)2317 PetscErrorCode PetscGaussLobattoLegendreElementLaplacianDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2318 {
2319 PetscErrorCode ierr;
2320
2321 PetscFunctionBegin;
2322 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr);
2323 ierr = PetscFree(*AA);CHKERRQ(ierr);
2324 *AA = NULL;
2325 PetscFunctionReturn(0);
2326 }
2327
2328 /*@C
2329 PetscGaussLobattoLegendreElementGradientCreate - computes the gradient for a single 1d GLL element
2330
2331 Not Collective
2332
2333 Input Parameter:
2334 + n - the number of GLL nodes
2335 . nodes - the GLL nodes
2336 . weights - the GLL weights
2337
2338 Output Parameter:
2339 . AA - the stiffness element
2340 - AAT - the transpose of AA (pass in NULL if you do not need this array)
2341
2342 Level: beginner
2343
2344 Notes:
2345 Destroy this with PetscGaussLobattoLegendreElementGradientDestroy()
2346
2347 You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented
2348
2349 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianDestroy()
2350
2351 @*/
PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal * nodes,PetscReal * weights,PetscReal *** AA,PetscReal *** AAT)2352 PetscErrorCode PetscGaussLobattoLegendreElementGradientCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
2353 {
2354 PetscReal **A, **AT = NULL;
2355 PetscErrorCode ierr;
2356 const PetscReal *gllnodes = nodes;
2357 const PetscInt p = n-1;
2358 PetscReal Li, Lj,d0;
2359 PetscInt i,j;
2360
2361 PetscFunctionBegin;
2362 ierr = PetscMalloc1(n,&A);CHKERRQ(ierr);
2363 ierr = PetscMalloc1(n*n,&A[0]);CHKERRQ(ierr);
2364 for (i=1; i<n; i++) A[i] = A[i-1]+n;
2365
2366 if (AAT) {
2367 ierr = PetscMalloc1(n,&AT);CHKERRQ(ierr);
2368 ierr = PetscMalloc1(n*n,&AT[0]);CHKERRQ(ierr);
2369 for (i=1; i<n; i++) AT[i] = AT[i-1]+n;
2370 }
2371
2372 if (n==1) {A[0][0] = 0.;}
2373 d0 = (PetscReal)p*((PetscReal)p+1.)/4.;
2374 for (i=0; i<n; i++) {
2375 for (j=0; j<n; j++) {
2376 A[i][j] = 0.;
2377 ierr = PetscDTComputeJacobi(0., 0., p, gllnodes[i], &Li);CHKERRQ(ierr);
2378 ierr = PetscDTComputeJacobi(0., 0., p, gllnodes[j], &Lj);CHKERRQ(ierr);
2379 if (i!=j) A[i][j] = Li/(Lj*(gllnodes[i]-gllnodes[j]));
2380 if ((j==i) && (i==0)) A[i][j] = -d0;
2381 if (j==i && i==p) A[i][j] = d0;
2382 if (AT) AT[j][i] = A[i][j];
2383 }
2384 }
2385 if (AAT) *AAT = AT;
2386 *AA = A;
2387 PetscFunctionReturn(0);
2388 }
2389
2390 /*@C
2391 PetscGaussLobattoLegendreElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGaussLobattoLegendreElementGradientCreate()
2392
2393 Not Collective
2394
2395 Input Parameter:
2396 + n - the number of GLL nodes
2397 . nodes - the GLL nodes
2398 . weights - the GLL weights
2399 . AA - the stiffness element
2400 - AAT - the transpose of the element
2401
2402 Level: beginner
2403
2404 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionCreate()
2405
2406 @*/
PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal * nodes,PetscReal * weights,PetscReal *** AA,PetscReal *** AAT)2407 PetscErrorCode PetscGaussLobattoLegendreElementGradientDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA,PetscReal ***AAT)
2408 {
2409 PetscErrorCode ierr;
2410
2411 PetscFunctionBegin;
2412 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr);
2413 ierr = PetscFree(*AA);CHKERRQ(ierr);
2414 *AA = NULL;
2415 if (*AAT) {
2416 ierr = PetscFree((*AAT)[0]);CHKERRQ(ierr);
2417 ierr = PetscFree(*AAT);CHKERRQ(ierr);
2418 *AAT = NULL;
2419 }
2420 PetscFunctionReturn(0);
2421 }
2422
2423 /*@C
2424 PetscGaussLobattoLegendreElementAdvectionCreate - computes the advection operator for a single 1d GLL element
2425
2426 Not Collective
2427
2428 Input Parameter:
2429 + n - the number of GLL nodes
2430 . nodes - the GLL nodes
2431 - weights - the GLL weightss
2432
2433 Output Parameter:
2434 . AA - the stiffness element
2435
2436 Level: beginner
2437
2438 Notes:
2439 Destroy this with PetscGaussLobattoLegendreElementAdvectionDestroy()
2440
2441 This is the same as the Gradient operator multiplied by the diagonal mass matrix
2442
2443 You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented
2444
2445 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementLaplacianCreate(), PetscGaussLobattoLegendreElementAdvectionDestroy()
2446
2447 @*/
PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal * nodes,PetscReal * weights,PetscReal *** AA)2448 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2449 {
2450 PetscReal **D;
2451 PetscErrorCode ierr;
2452 const PetscReal *gllweights = weights;
2453 const PetscInt glln = n;
2454 PetscInt i,j;
2455
2456 PetscFunctionBegin;
2457 ierr = PetscGaussLobattoLegendreElementGradientCreate(n,nodes,weights,&D,NULL);CHKERRQ(ierr);
2458 for (i=0; i<glln; i++){
2459 for (j=0; j<glln; j++) {
2460 D[i][j] = gllweights[i]*D[i][j];
2461 }
2462 }
2463 *AA = D;
2464 PetscFunctionReturn(0);
2465 }
2466
2467 /*@C
2468 PetscGaussLobattoLegendreElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element
2469
2470 Not Collective
2471
2472 Input Parameter:
2473 + n - the number of GLL nodes
2474 . nodes - the GLL nodes
2475 . weights - the GLL weights
2476 - A - advection
2477
2478 Level: beginner
2479
2480 .seealso: PetscDTGaussLobattoLegendreQuadrature(), PetscGaussLobattoLegendreElementAdvectionCreate()
2481
2482 @*/
PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal * nodes,PetscReal * weights,PetscReal *** AA)2483 PetscErrorCode PetscGaussLobattoLegendreElementAdvectionDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2484 {
2485 PetscErrorCode ierr;
2486
2487 PetscFunctionBegin;
2488 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr);
2489 ierr = PetscFree(*AA);CHKERRQ(ierr);
2490 *AA = NULL;
2491 PetscFunctionReturn(0);
2492 }
2493
PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal * nodes,PetscReal * weights,PetscReal *** AA)2494 PetscErrorCode PetscGaussLobattoLegendreElementMassCreate(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2495 {
2496 PetscReal **A;
2497 PetscErrorCode ierr;
2498 const PetscReal *gllweights = weights;
2499 const PetscInt glln = n;
2500 PetscInt i,j;
2501
2502 PetscFunctionBegin;
2503 ierr = PetscMalloc1(glln,&A);CHKERRQ(ierr);
2504 ierr = PetscMalloc1(glln*glln,&A[0]);CHKERRQ(ierr);
2505 for (i=1; i<glln; i++) A[i] = A[i-1]+glln;
2506 if (glln==1) {A[0][0] = 0.;}
2507 for (i=0; i<glln; i++) {
2508 for (j=0; j<glln; j++) {
2509 A[i][j] = 0.;
2510 if (j==i) A[i][j] = gllweights[i];
2511 }
2512 }
2513 *AA = A;
2514 PetscFunctionReturn(0);
2515 }
2516
PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal * nodes,PetscReal * weights,PetscReal *** AA)2517 PetscErrorCode PetscGaussLobattoLegendreElementMassDestroy(PetscInt n,PetscReal *nodes,PetscReal *weights,PetscReal ***AA)
2518 {
2519 PetscErrorCode ierr;
2520
2521 PetscFunctionBegin;
2522 ierr = PetscFree((*AA)[0]);CHKERRQ(ierr);
2523 ierr = PetscFree(*AA);CHKERRQ(ierr);
2524 *AA = NULL;
2525 PetscFunctionReturn(0);
2526 }
2527
2528 /*@
2529 PetscDTIndexToBary - convert an index into a barycentric coordinate.
2530
2531 Input Parameters:
2532 + 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)
2533 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
2534 - index - the index to convert: should be >= 0 and < Binomial(len - 1 + sum, sum)
2535
2536 Output Parameter:
2537 . coord - will be filled with the barycentric coordinate
2538
2539 Level: beginner
2540
2541 Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the
2542 least significant and the last index is the most significant.
2543
2544 .seealso: PetscDTBaryToIndex()
2545 @*/
PetscDTIndexToBary(PetscInt len,PetscInt sum,PetscInt index,PetscInt coord[])2546 PetscErrorCode PetscDTIndexToBary(PetscInt len, PetscInt sum, PetscInt index, PetscInt coord[])
2547 {
2548 PetscInt c, d, s, total, subtotal, nexttotal;
2549
2550 PetscFunctionBeginHot;
2551 if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
2552 if (index < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index must be non-negative");
2553 if (!len) {
2554 if (!sum && !index) PetscFunctionReturn(0);
2555 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
2556 }
2557 for (c = 1, total = 1; c <= len; c++) {
2558 /* total is the number of ways to have a tuple of length c with sum */
2559 if (index < total) break;
2560 total = (total * (sum + c)) / c;
2561 }
2562 if (c > len) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index out of range");
2563 for (d = c; d < len; d++) coord[d] = 0;
2564 for (s = 0, subtotal = 1, nexttotal = 1; c > 0;) {
2565 /* subtotal is the number of ways to have a tuple of length c with sum s */
2566 /* nexttotal is the number of ways to have a tuple of length c-1 with sum s */
2567 if ((index + subtotal) >= total) {
2568 coord[--c] = sum - s;
2569 index -= (total - subtotal);
2570 sum = s;
2571 total = nexttotal;
2572 subtotal = 1;
2573 nexttotal = 1;
2574 s = 0;
2575 } else {
2576 subtotal = (subtotal * (c + s)) / (s + 1);
2577 nexttotal = (nexttotal * (c - 1 + s)) / (s + 1);
2578 s++;
2579 }
2580 }
2581 PetscFunctionReturn(0);
2582 }
2583
2584 /*@
2585 PetscDTBaryToIndex - convert a barycentric coordinate to an index
2586
2587 Input Parameters:
2588 + 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)
2589 . sum - the value that the sum of the barycentric coordinates (which will be non-negative integers) should sum to
2590 - coord - a barycentric coordinate with the given length and sum
2591
2592 Output Parameter:
2593 . index - the unique index for the coordinate, >= 0 and < Binomial(len - 1 + sum, sum)
2594
2595 Level: beginner
2596
2597 Note: the indices map to barycentric coordinates in lexicographic order, where the first index is the
2598 least significant and the last index is the most significant.
2599
2600 .seealso: PetscDTIndexToBary
2601 @*/
PetscDTBaryToIndex(PetscInt len,PetscInt sum,const PetscInt coord[],PetscInt * index)2602 PetscErrorCode PetscDTBaryToIndex(PetscInt len, PetscInt sum, const PetscInt coord[], PetscInt *index)
2603 {
2604 PetscInt c;
2605 PetscInt i;
2606 PetscInt total;
2607
2608 PetscFunctionBeginHot;
2609 if (len < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length must be non-negative");
2610 if (!len) {
2611 if (!sum) {
2612 *index = 0;
2613 PetscFunctionReturn(0);
2614 }
2615 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid index or sum for length 0 barycentric coordinate");
2616 }
2617 for (c = 1, total = 1; c < len; c++) total = (total * (sum + c)) / c;
2618 i = total - 1;
2619 c = len - 1;
2620 sum -= coord[c];
2621 while (sum > 0) {
2622 PetscInt subtotal;
2623 PetscInt s;
2624
2625 for (s = 1, subtotal = 1; s < sum; s++) subtotal = (subtotal * (c + s)) / s;
2626 i -= subtotal;
2627 sum -= coord[--c];
2628 }
2629 *index = i;
2630 PetscFunctionReturn(0);
2631 }
2632