Intrepid2
Intrepid2_PolylibDef.hpp
Go to the documentation of this file.
1/*
2// @HEADER
3// *****************************************************************************
4// Intrepid2 Package
5//
6// Copyright 2007 NTESS and the Intrepid2 contributors.
7// SPDX-License-Identifier: BSD-3-Clause
8// *****************************************************************************
9// @HEADER
10*/
11
13//
14// File: Intrepid_PolylibDef.hpp
15//
16// For more information, please see: http://www.nektar.info
17//
18// The MIT License
19//
20// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
21// Department of Aeronautics, Imperial College London (UK), and Scientific
22// Computing and Imaging Institute, University of Utah (USA).
23//
24// License for the specific language governing rights and limitations under
25// Permission is hereby granted, free of charge, to any person obtaining a
26// copy of this software and associated documentation files (the "Software"),
27// to deal in the Software without restriction, including without limitation
28// the rights to use, copy, modify, merge, publish, distribute, sublicense,
29// and/or sell copies of the Software, and to permit persons to whom the
30// Software is furnished to do so, subject to the following conditions:
31//
32// The above copyright notice and this permission notice shall be included
33// in all copies or substantial portions of the Software.
34//
35// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
36// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
37// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
38// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
39// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
40// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
41// DEALINGS IN THE SOFTWARE.
42//
43// Description:
44// This file is redistributed with the Intrepid package. It should be used
45// in accordance with the above MIT license, at the request of the authors.
46// This file is NOT covered by the usual Intrepid/Trilinos LGPL license.
47//
48// Origin: Nektar++ library, http://www.nektar.info, downloaded on
49// March 10, 2009.
50//
52
53
62#ifndef __INTREPID2_POLYLIB_DEF_HPP__
63#define __INTREPID2_POLYLIB_DEF_HPP__
64
65#if defined(_MSC_VER) || defined(_WIN32) && defined(__ICL)
66// M_PI, M_SQRT2, etc. are hidden in MSVC by #ifdef _USE_MATH_DEFINES
67 #ifndef _USE_MATH_DEFINES
68 #define _USE_MATH_DEFINES
69 #endif
70 #include <math.h>
71#endif
72
73namespace Intrepid2 {
74
75 // -----------------------------------------------------------------------
76 // Points and Weights
77 // -----------------------------------------------------------------------
78
79 template<>
80 template<typename zViewType,
81 typename wViewType>
82 KOKKOS_INLINE_FUNCTION
83 void
84 Polylib::Serial::Cubature<POLYTYPE_GAUSS>::
85 getValues( zViewType z,
86 wViewType w,
87 const ordinal_type np,
88 const double alpha,
89 const double beta) {
90 const double one = 1.0, two = 2.0, apb = alpha + beta;
91 double fac;
92
93 JacobiZeros(z, np, alpha, beta);
94 JacobiPolynomialDerivative(np, z, w, np, alpha, beta);
95
96 fac = pow(two, apb + one)*GammaFunction(alpha + np + one)*GammaFunction(beta + np + one);
97 fac /= GammaFunction((np + one))*GammaFunction(apb + np + one);
98
99 for (ordinal_type i = 0; i < np; ++i)
100 w(i) = fac/(w(i)*w(i)*(one-z(i)*z(i)));
101 }
102
103 template<>
104 template<typename zViewType,
105 typename wViewType>
106 KOKKOS_INLINE_FUNCTION
107 void
108 Polylib::Serial::Cubature<POLYTYPE_GAUSS_RADAU_LEFT>::
109 getValues( zViewType z,
110 wViewType w,
111 const ordinal_type np,
112 const double alpha,
113 const double beta) {
114 if (np == 1) {
115 z(0) = 0.0;
116 w(0) = 2.0;
117 } else {
118 const double one = 1.0, two = 2.0, apb = alpha + beta;
119 double fac;
120
121 z(0) = -one;
122
123 auto z_plus_1 = Kokkos::subview(z, Kokkos::pair<ordinal_type,ordinal_type>(1, z.extent(0)));
124 JacobiZeros(z_plus_1, np-1, alpha, beta+1);
125
126 Kokkos::View<typename zViewType::value_type*,Kokkos::AnonymousSpace,Kokkos::MemoryUnmanaged> null;
127 JacobiPolynomial(np, z, w, null, np-1, alpha, beta);
128
129 fac = pow(two, apb)*GammaFunction(alpha + np)*GammaFunction(beta + np);
130 fac /= GammaFunction(np)*(beta + np)*GammaFunction(apb + np + 1);
131
132 for (ordinal_type i = 0; i < np; ++i)
133 w(i) = fac*(1-z(i))/(w(i)*w(i));
134 w(0) *= (beta + one);
135 }
136 }
137
138 template<>
139 template<typename zViewType,
140 typename wViewType>
141 KOKKOS_INLINE_FUNCTION
142 void
143 Polylib::Serial::Cubature<POLYTYPE_GAUSS_RADAU_RIGHT>::
144 getValues( zViewType z,
145 wViewType w,
146 const ordinal_type np,
147 const double alpha,
148 const double beta) {
149 if (np == 1) {
150 z(0) = 0.0;
151 w(0) = 2.0;
152 } else {
153 const double one = 1.0, two = 2.0, apb = alpha + beta;
154 double fac;
155
156 JacobiZeros(z, np-1, alpha+1, beta);
157 z(np-1) = one;
158
159 Kokkos::View<typename zViewType::value_type*,Kokkos::AnonymousSpace,Kokkos::MemoryUnmanaged> null;
160 JacobiPolynomial(np, z, w, null, np-1, alpha, beta);
161
162 fac = pow(two,apb)*GammaFunction(alpha + np)*GammaFunction(beta + np);
163 fac /= GammaFunction(np)*(alpha + np)*GammaFunction(apb + np + 1);
164
165 for (ordinal_type i = 0; i < np; ++i)
166 w(i) = fac*(1+z(i))/(w(i)*w(i));
167 w(np-1) *= (alpha + one);
168 }
169 }
170
171 template<>
172 template<typename zViewType,
173 typename wViewType>
174 KOKKOS_INLINE_FUNCTION
175 void
176 Polylib::Serial::Cubature<POLYTYPE_GAUSS_LOBATTO>::
177 getValues( zViewType z,
178 wViewType w,
179 const ordinal_type np,
180 const double alpha,
181 const double beta) {
182 if ( np == 1 ) {
183 z(0) = 0.0;
184 w(0) = 2.0;
185 } else {
186 const double one = 1.0, apb = alpha + beta, two = 2.0;
187 double fac;
188
189 z(0) = -one;
190 z(np-1) = one;
191
192 auto z_plus_1 = Kokkos::subview(z, Kokkos::pair<ordinal_type,ordinal_type>(1, z.extent(0)));
193 JacobiZeros(z_plus_1, np-2, alpha+one, beta+one);
194
195 Kokkos::View<typename zViewType::value_type*,Kokkos::AnonymousSpace,Kokkos::MemoryUnmanaged> null;
196 JacobiPolynomial(np, z, w, null, np-1, alpha, beta);
197
198 fac = pow(two, apb + 1)*GammaFunction(alpha + np)*GammaFunction(beta + np);
199 fac /= (np-1)*GammaFunction(np)*GammaFunction(alpha + beta + np + one);
200
201 for (ordinal_type i = 0; i < np; ++i)
202 w(i) = fac/(w(i)*w(i));
203
204 w(0) *= (beta + one);
205 w(np-1) *= (alpha + one);
206 }
207 }
208
209 // -----------------------------------------------------------------------
210 // Derivatives
211 // -----------------------------------------------------------------------
212
213 template<>
214 template<typename DViewType,
215 typename zViewType>
216 KOKKOS_INLINE_FUNCTION
217 void
218 Polylib::Serial::Derivative<POLYTYPE_GAUSS>::
219 getValues( DViewType D,
220 const zViewType z,
221 const ordinal_type np,
222 const double alpha,
223 const double beta) {
224 if (np <= 0) {
225 D(0,0) = 0.0;
226 } else {
227 const double one = 1.0, two = 2.0;
228
229 auto pd = Kokkos::subview(D, np-1, Kokkos::pair<int,int>(0,np));
230 JacobiPolynomialDerivative(np, z, pd, np, alpha, beta);
231
232 // The temporary view pd is stored in the last row of the matrix D
233 // This loop is designed so that we do not overwrite pd entries before we read them
234 for (ordinal_type i = 0; i < np; ++i) {
235 const auto & pd_i = pd(i);
236 const auto & z_i = z(i);
237 for (ordinal_type j = 0; j < i; ++j) {
238 const auto & pd_j = pd(j);
239 const auto & z_j = z(j);
240 D(j,i) = pd_j/(pd_i*(z_j-z_i));
241 D(i,j) = pd_i/(pd_j*(z_i-z_j));
242 }
243 D(i,i) = (alpha - beta + (alpha + beta + two)*z_i) / (two*(one - z_i*z_i));
244 }
245 }
246 }
247
248 template<>
249 template<typename DViewType,
250 typename zViewType>
251 KOKKOS_INLINE_FUNCTION
252 void
253 Polylib::Serial::Derivative<POLYTYPE_GAUSS_RADAU_LEFT>::
254 getValues( DViewType D,
255 const zViewType z,
256 const ordinal_type np,
257 const double alpha,
258 const double beta) {
259 if (np <= 0) {
260 D(0,0) = 0.0;
261 } else {
262 const double one = 1.0, two = 2.0;
263
264 auto pd = Kokkos::subview(D, np-1, Kokkos::pair<int,int>(0,np));
265 pd(0) = pow(-one,np-1)*GammaFunction(np+beta+one) / (GammaFunction(np)*GammaFunction(beta+two));
266
267 auto pd_plus_1 = Kokkos::subview(pd, Kokkos::pair<ordinal_type,ordinal_type>(1, pd.extent(0)));
268 auto z_plus_1 = Kokkos::subview( z, Kokkos::pair<ordinal_type,ordinal_type>(1, z.extent(0)));
269
270 JacobiPolynomialDerivative(np-1, z_plus_1, pd_plus_1, np-1, alpha, beta+1);
271 for(ordinal_type i = 1; i < np; ++i)
272 pd(i) *= (1+z(i));
273
274 // The temporary view pd is stored in the last row of the matrix D
275 // This loop is designed so that we do not overwrite pd entries before we read them
276 for (ordinal_type i = 0; i < np; ++i) {
277 const auto & pd_i = pd(i);
278 const auto & z_i = z(i);
279 for (ordinal_type j = 0; j < i; ++j) {
280 const auto & pd_j = pd(j);
281 const auto & z_j = z(j);
282 D(j,i) = pd_j/(pd_i*(z_j-z_i));
283 D(i,j) = pd_i/(pd_j*(z_i-z_j));
284 }
285 if (i == 0)
286 D(i,i) = -(np + alpha + beta + one)*(np - one) / (two*(beta + two));
287 else
288 D(i,i) = (alpha - beta + one + (alpha + beta + one)*z_i) / (two*(one - z_i*z_i));
289 }
290 }
291 }
292
293 template<>
294 template<typename DViewType,
295 typename zViewType>
296 KOKKOS_INLINE_FUNCTION
297 void
298 Polylib::Serial::Derivative<POLYTYPE_GAUSS_RADAU_RIGHT>::
299 getValues( DViewType D,
300 const zViewType z,
301 const ordinal_type np,
302 const double alpha,
303 const double beta) {
304 if (np <= 0) {
305 D(0,0) = 0.0;
306 } else {
307 const double one = 1.0, two = 2.0;
308
309 auto pd = Kokkos::subview(D, np-1, Kokkos::pair<int,int>(0,np));
310
311 JacobiPolynomialDerivative(np-1, z, pd, np-1, alpha+1, beta);
312 for (ordinal_type i = 0; i < np-1; ++i)
313 pd(i) *= (1-z(i));
314
315 pd(np-1) = -GammaFunction(np+alpha+one) / (GammaFunction(np)*GammaFunction(alpha+two));
316
317 // The temporary view pd is stored in the last row of the matrix D
318 // This loop is designed so that we do not overwrite pd entries before we read them
319 for (ordinal_type i = 0; i < np; ++i) {
320 const auto & pd_i = pd(i);
321 const auto & z_i = z(i);
322 for (ordinal_type j = 0; j < i; ++j) {
323 const auto & pd_j = pd(j);
324 const auto & z_j = z(j);
325 D(j,i) = pd_j/(pd_i*(z_j-z_i));
326 D(i,j) = pd_i/(pd_j*(z_i-z_j));
327 }
328 if (i == np-1)
329 D(i,i) = (np + alpha + beta + one)*(np - one) / (two*(alpha + two));
330 else
331 D(i,i) = (alpha - beta - one + (alpha + beta + one)*z_i) / (two*(one - z_i*z_i));
332 }
333 }
334 }
335
336 template<>
337 template<typename DViewType,
338 typename zViewType>
339 KOKKOS_INLINE_FUNCTION
340 void
341 Polylib::Serial::Derivative<POLYTYPE_GAUSS_LOBATTO>::
342 getValues( DViewType D,
343 const zViewType z,
344 const ordinal_type np,
345 const double alpha,
346 const double beta) {
347 if (np <= 0) {
348 D(0,0) = 0.0;
349 } else {
350 const double one = 1.0, two = 2.0;
351
352 auto pd = Kokkos::subview(D, np-1, Kokkos::pair<int,int>(0,np));
353
354 pd(0) = two*pow(-one,np)*GammaFunction(np + beta);
355 pd(0) /= GammaFunction(np - one)*GammaFunction(beta + two);
356
357 auto pd_plus_1 = Kokkos::subview(pd, Kokkos::pair<ordinal_type,ordinal_type>(1, pd.extent(0)));
358 auto z_plus_1 = Kokkos::subview( z, Kokkos::pair<ordinal_type,ordinal_type>(1, z.extent(0)));
359
360 JacobiPolynomialDerivative(np-2, z_plus_1, pd_plus_1, np-2, alpha+1, beta+1);
361 for (ordinal_type i = 1; i < np-1; ++i) {
362 const auto & z_i = z(i);
363 pd(i) *= (one-z_i*z_i);
364 }
365
366 pd(np-1) = -two*GammaFunction(np + alpha);
367 pd(np-1) /= GammaFunction(np - one)*GammaFunction(alpha + two);
368
369 // The temporary view pd is stored in the last row of the matrix D
370 // This loop is designed so that we do not overwrite pd entries before we read them
371 for (ordinal_type i = 0; i < np; ++i) {
372 const auto & pd_i = pd(i);
373 const auto & z_i = z(i);
374 for (ordinal_type j = 0; j < i; ++j) {
375 const auto & pd_j = pd(j);
376 const auto & z_j = z(j);
377 D(j,i) = pd_j/(pd_i*(z_j-z_i));
378 D(i,j) = pd_i/(pd_j*(z_i-z_j));
379 }
380 if (i == 0)
381 D(i,i) = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
382 else if (i == np-1)
383 D(i,i) =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
384 else
385 D(i,i) = (alpha - beta + (alpha + beta)*z_i)/(two*(one - z_i*z_i));
386 }
387 }
388 }
389
390 // -----------------------------------------------------------------------
391 // Lagrangian Interpolants
392 // -----------------------------------------------------------------------
393
394 template<>
395 template<typename zViewType>
396 KOKKOS_INLINE_FUNCTION
397 typename zViewType::value_type
398 Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS>::
399 getValue(const ordinal_type i,
400 const typename zViewType::value_type z,
401 const zViewType zgj,
402 const ordinal_type np,
403 const double alpha,
404 const double beta) {
405 const double tol = tolerence();
406
407 typedef typename zViewType::value_type value_type;
408 typedef typename zViewType::pointer_type pointer_type;
409
410 value_type h, p_buf, pd_buf, zi_buf = zgj(i);
411 Kokkos::View<value_type*,Kokkos::AnonymousSpace,Kokkos::MemoryUnmanaged>
412 p((pointer_type)&p_buf, 1), pd((pointer_type)&pd_buf, 1), zi((pointer_type)&zi_buf, 1),
413 zv(const_cast<pointer_type>(&z), 1), null;
414
415 const auto dz = z - zi(0);
416 if (Util<value_type>::abs(dz) < tol)
417 return 1.0;
418
419 JacobiPolynomialDerivative(1, zi, pd, np, alpha, beta);
420 JacobiPolynomial(1, zv, p, null , np, alpha, beta);
421
422 h = p(0)/(pd(0)*dz);
423
424 return h;
425 }
426
427 template<>
428 template<typename zViewType>
429 KOKKOS_INLINE_FUNCTION
430 typename zViewType::value_type
431 Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_RADAU_LEFT>::
432 getValue(const ordinal_type i,
433 const typename zViewType::value_type z,
434 const zViewType zgrj,
435 const ordinal_type np,
436 const double alpha,
437 const double beta) {
438 const double tol = tolerence();
439
440 typedef typename zViewType::value_type value_type;
441 typedef typename zViewType::pointer_type pointer_type;
442
443 value_type h, p_buf, pd_buf, zi_buf = zgrj(i);
444 Kokkos::View<value_type*,Kokkos::AnonymousSpace,Kokkos::MemoryUnmanaged>
445 p((pointer_type)&p_buf, 1), pd((pointer_type)&pd_buf, 1), zi((pointer_type)&zi_buf, 1),
446 zv(const_cast<pointer_type>(&z), 1), null;
447
448 const auto dz = z - zi(0);
449 if (Util<value_type>::abs(dz) < tol)
450 return 1.0;
451
452 JacobiPolynomial(1, zi, p , null, np-1, alpha, beta + 1);
453
454 // need to use this routine in case zi = -1 or 1
455 JacobiPolynomialDerivative(1, zi, pd, np-1, alpha, beta + 1);
456 h = (1.0 + zi(0))*pd(0) + p(0);
457
458 JacobiPolynomial(1, zv, p, null, np-1, alpha, beta + 1);
459 h = (1.0 + z)*p(0)/(h*dz);
460
461 return h;
462 }
463
464
465 template<>
466 template<typename zViewType>
467 KOKKOS_INLINE_FUNCTION
468 typename zViewType::value_type
469 Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_RADAU_RIGHT>::
470 getValue(const ordinal_type i,
471 const typename zViewType::value_type z,
472 const zViewType zgrj,
473 const ordinal_type np,
474 const double alpha,
475 const double beta) {
476 const double tol = tolerence();
477
478 typedef typename zViewType::value_type value_type;
479 typedef typename zViewType::pointer_type pointer_type;
480
481 value_type h, p_buf, pd_buf, zi_buf = zgrj(i);
482 Kokkos::View<value_type*,Kokkos::AnonymousSpace,Kokkos::MemoryUnmanaged>
483 p((pointer_type)&p_buf, 1), pd((pointer_type)&pd_buf, 1), zi((pointer_type)&zi_buf, 1),
484 zv(const_cast<pointer_type>(&z), 1), null;
485
486 const auto dz = z - zi(0);
487 if (Util<value_type>::abs(dz) < tol)
488 return 1.0;
489
490 JacobiPolynomial(1, zi, p , null, np-1, alpha+1, beta);
491
492 // need to use this routine in case z = -1 or 1
493 JacobiPolynomialDerivative(1, zi, pd, np-1, alpha+1, beta);
494 h = (1.0 - zi(0))*pd(0) - p(0);
495
496 JacobiPolynomial (1, zv, p, null, np-1, alpha+1, beta);
497 h = (1.0 - z)*p(0)/(h*dz);
498
499 return h;
500 }
501
502
503 template<>
504 template<typename zViewType>
505 KOKKOS_INLINE_FUNCTION
506 typename zViewType::value_type
507 Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_LOBATTO>::
508 getValue(const ordinal_type i,
509 const typename zViewType::value_type z,
510 const zViewType zglj,
511 const ordinal_type np,
512 const double alpha,
513 const double beta) {
514 const double one = 1.0, two = 2.0, tol = tolerence();
515
516 typedef typename zViewType::value_type value_type;
517 typedef typename zViewType::pointer_type pointer_type;
518
519 value_type h, p_buf, pd_buf, zi_buf = zglj(i);
520 Kokkos::View<value_type*,Kokkos::AnonymousSpace,Kokkos::MemoryUnmanaged>
521 p((pointer_type)&p_buf, 1), pd((pointer_type)&pd_buf, 1), zi((pointer_type)&zi_buf, 1),
522 zv(const_cast<pointer_type>(&z), 1), null;
523
524 const auto dz = z - zi(0);
525 if (Util<value_type>::abs(dz) < tol)
526 return 1.0;
527
528 JacobiPolynomial(1, zi, p , null, np-2, alpha + one, beta + one);
529
530 // need to use this routine in case z = -1 or 1
531 JacobiPolynomialDerivative(1, zi, pd, np-2, alpha + one, beta + one);
532 h = (one - zi(0)*zi(0))*pd(0) - two*zi(0)*p(0);
533
534 JacobiPolynomial(1, zv, p, null, np-2, alpha + one, beta + one);
535 h = (one - z*z)*p(0)/(h*dz);
536
537 return h;
538 }
539
540 // -----------------------------------------------------------------------
541 // Interpolation Operator
542 // -----------------------------------------------------------------------
543
544 template<EPolyType polyType>
545 template<typename imViewType,
546 typename zgrjViewType,
547 typename zmViewType>
548 KOKKOS_INLINE_FUNCTION
549 void
550 Polylib::Serial::InterpolationOperator<polyType>::
551 getMatrix( imViewType im,
552 const zgrjViewType zgrj,
553 const zmViewType zm,
554 const ordinal_type nz,
555 const ordinal_type mz,
556 const double alpha,
557 const double beta) {
558 for (ordinal_type i = 0; i < mz; ++i) {
559 const auto zp = zm(i);
560 for (ordinal_type j = 0; j < nz; ++j)
561 im(i, j) = LagrangianInterpolant<polyType>::getValue(j, zp, zgrj, nz, alpha, beta);
562 }
563 }
564
565 // -----------------------------------------------------------------------
566
567 template<typename zViewType,
568 typename polyiViewType,
569 typename polydViewType>
570 KOKKOS_INLINE_FUNCTION
571 void
572 Polylib::Serial::
573 JacobiPolynomial(const ordinal_type np,
574 const zViewType z,
575 polyiViewType polyi,
576 polydViewType polyd,
577 const ordinal_type n,
578 const double alpha,
579 const double beta) {
580 const double zero = 0.0, one = 1.0, two = 2.0;
581
582 if (! np) {
583 return;
584 }
585
586 if (n == 0) {
587 if (polyi.data())
588 for (ordinal_type i = 0; i < np; ++i)
589 polyi(i) = one;
590 if (polyd.data())
591 for (ordinal_type i = 0; i < np; ++i)
592 polyd(i) = zero;
593 } else if (n == 1) {
594 if (polyi.data())
595 for (ordinal_type i = 0; i < np; ++i)
596 polyi(i) = 0.5*(alpha - beta + (alpha + beta + two)*z(i));
597 if (polyd.data())
598 for (ordinal_type i = 0; i < np; ++i)
599 polyd(i) = 0.5*(alpha + beta + two);
600 } else {
601 INTREPID2_TEST_FOR_ABORT(polyd.data() && !polyd.data() ,
602 ">>> ERROR (Polylib::Serial::JacobiPolynomial): polyi view needed to compute polyd view.");
603 if(!polyi.data()) return;
604
605 constexpr ordinal_type maxOrder = 2*MaxPolylibPoint-1;
606
607 INTREPID2_TEST_FOR_ABORT(maxOrder < n,
608 ">>> ERROR (Polylib::Serial::JacobiPolynomial): Requested order exceeds maxOrder .");
609
610 double a2[maxOrder-1]={}, a3[maxOrder-1]={}, a4[maxOrder-1]={};
611 double ad1(0.0), ad2(0.0), ad3(0.0);
612 const double apb = alpha + beta;
613 const double amb = alpha - beta;
614
615
616 for (auto k = 2; k <= n; ++k) {
617 double a1 = two*k*(k + apb)*(two*k + apb - two);
618 a2[k-2] = (two*k + apb - one)*(apb*amb)/a1;
619 a3[k-2] = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb)/a1;
620 a4[k-2] = two*(k + alpha - one)*(k + beta - one)*(two*k + apb)/a1;
621 }
622
623 if (polyd.data()) {
624 double ad4 = (two*n + alpha + beta);
625 ad1 = n*(alpha - beta)/ad4;
626 ad2 = n*(two*n + alpha + beta)/ad4;
627 ad3 = two*(n + alpha)*(n + beta)/ad4;
628 }
629
630 typename polyiViewType::value_type polyn0, polyn1, polyn2;
631 for(ordinal_type i = 0; i < np; ++i) {
632 const auto & z_i = z(i);
633 polyn2 = one;
634 polyn1 = 0.5*(amb + (apb + two)*z_i);
635 polyn0 = (a2[0] + a3[0]*z_i)*polyn1 - a4[0]*polyn2;
636 for (ordinal_type k = 1; k < n-1; ++k) {
637 polyn2 = polyn1;
638 polyn1 = polyn0;
639 polyn0 = (a2[k] + a3[k]*z_i)*polyn1 - a4[k]*polyn2;
640 }
641 if (polyd.data()) {
642 polyd(i) = (ad1- ad2*z_i)*polyn0 + ad3*polyn1 / (one - z_i*z_i);
643 }
644 polyi(i) = polyn0;
645 }
646 }
647 }
648
649 template<typename zViewType,
650 typename polydViewType>
651 KOKKOS_INLINE_FUNCTION
652 void
653 Polylib::Serial::
654 JacobiPolynomialDerivative(const ordinal_type np,
655 const zViewType z,
656 polydViewType polyd,
657 const ordinal_type n,
658 const double alpha,
659 const double beta) {
660 const double one = 1.0;
661 if (n == 0)
662 for(ordinal_type i = 0; i < np; ++i)
663 polyd(i) = 0.0;
664 else {
665 Kokkos::View<typename polydViewType::value_type*,Kokkos::AnonymousSpace,Kokkos::MemoryUnmanaged> null;
666 JacobiPolynomial(np, z, polyd, null, n-1, alpha+one, beta+one);
667 for(ordinal_type i = 0; i < np; ++i)
668 polyd(i) *= 0.5*(alpha + beta + n + one);
669 }
670 }
671
672 // -----------------------------------------------------------------------
673
674 template<typename zViewType,
675 bool DeflationEnabled>
676 KOKKOS_INLINE_FUNCTION
677 void
678 Polylib::Serial::
679 JacobiZeros( zViewType z,
680 const ordinal_type n,
681 const double alpha,
682 const double beta) {
683 if (DeflationEnabled)
684 JacobiZerosPolyDeflation(z, n, alpha, beta);
685 else
686 JacobiZerosTriDiagonal(z, n, alpha, beta);
687 }
688
689 template<typename zViewType>
690 KOKKOS_INLINE_FUNCTION
691 void
692 Polylib::Serial::
693 JacobiZerosPolyDeflation( zViewType z,
694 const ordinal_type n,
695 const double alpha,
696 const double beta) {
697 if(!n)
698 return;
699
700 const double dth = M_PI/(2.0*n);
701 const double one = 1.0, two = 2.0;
702 const double tol = tolerence();
703
704 typedef typename zViewType::value_type value_type;
705 typedef typename zViewType::pointer_type pointer_type;
706
707 value_type r_buf, poly_buf, pder_buf;
708 Kokkos::View<value_type*,Kokkos::AnonymousSpace,Kokkos::MemoryUnmanaged>
709 poly((pointer_type)&poly_buf, 1), pder((pointer_type)&pder_buf, 1), r((pointer_type)&r_buf, 1);
710
711 value_type rlast = 0.0;
712 for (auto k = 0; k < n; ++k) {
713 r(0) = -cos((two*(double)k + one) * dth);
714 if (k)
715 r(0) = 0.5*(r(0) + rlast);
716
717 for (ordinal_type j = 1; j < MaxPolylibIteration; ++j) {
718 JacobiPolynomial(1, r, poly, pder, n, alpha, beta);
719
720 value_type sum = 0;
721 for (ordinal_type i = 0; i < k; ++i)
722 sum += one/(r(0) - z(i));
723
724 const value_type delr = -poly(0) / (pder(0) - sum * poly(0));
725 r(0) += delr;
726
727 if( Util<value_type>::abs(delr) < tol )
728 break;
729 }
730 z(k) = r(0);
731 rlast = r(0);
732 }
733 }
734
735 template<typename aViewType>
736 KOKKOS_INLINE_FUNCTION
737 void
738 Polylib::Serial::
739 JacobiZerosTriDiagonal( aViewType a,
740 const ordinal_type n,
741 const double alpha,
742 const double beta) {
743 if(!n)
744 return;
745
746 typedef typename aViewType::value_type value_type;
747 typedef typename aViewType::pointer_type pointer_type;
748
749 value_type b_buf[MaxPolylibPoint];
750 Kokkos::View<value_type*,Kokkos::AnonymousSpace,Kokkos::MemoryUnmanaged>
751 b((pointer_type)&b_buf[0], MaxPolylibPoint);
752
753 // generate normalised terms
754 const double apb = alpha + beta;
755 double apbi = 2.0 + apb;
756
757 b(n-1) = pow(2.0,apb+1.0)*GammaFunction(alpha+1.0)*GammaFunction(beta+1.0)/GammaFunction(apbi);
758 a(0) = (beta-alpha)/apbi;
759 b(0) = sqrt(4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
760
761 const double a2b2 = beta*beta-alpha*alpha;
762 for (ordinal_type i = 1; i < n-1; ++i) {
763 apbi = 2.0*(i+1) + apb;
764 a(i) = a2b2/((apbi-2.0)*apbi);
765 b(i) = sqrt(4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
766 ((apbi*apbi-1)*apbi*apbi));
767 }
768
769 apbi = 2.0*n + apb;
770 //a(n-1) = a2b2/((apbi-2.0)*apbi); // THIS IS A BUG!!!
771 if (n>1) a(n-1) = a2b2/((apbi-2.0)*apbi);
772
773 // find eigenvalues
774 TriQL(a, b, n);
775 }
776
777
778 template<typename dViewType,
779 typename eViewType>
780 KOKKOS_INLINE_FUNCTION
781 void
782 Polylib::Serial::
783 TriQL( dViewType d,
784 eViewType e,
785 const ordinal_type n) {
786 ordinal_type m,l,iter,i,k;
787
788 typedef typename dViewType::value_type value_type;
789 value_type s,r,p,g,f,dd,c,b;
790
791 for (l=0; l<n; ++l) {
792 iter=0;
793 do {
794 for (m=l; m<n-1; ++m) {
795 dd=Util<value_type>::abs(d(m))+Util<value_type>::abs(d(m+1));
796 if (Util<value_type>::abs(e(m))+dd == dd) break;
797 }
798 if (m != l) {
799 if (iter++ == MaxPolylibIteration) {
800 INTREPID2_TEST_FOR_ABORT(true,
801 ">>> ERROR (Polylib::Serial): Too many iterations in TQLI.");
802 }
803 g=(d(l+1)-d(l))/(2.0*e(l));
804 r=sqrt((g*g)+1.0);
805 //g=d(m)-d(l)+e(l)/(g+sign(r,g));
806 g=d(m)-d(l)+e(l)/(g+((g)<0 ? value_type(-Util<value_type>::abs(r)) : Util<value_type>::abs(r)));
807 s=c=1.0;
808 p=0.0;
809 for (i=m-1; i>=l; i--) {
810 f=s*e(i);
811 b=c*e(i);
812 if (Util<value_type>::abs(f) >= Util<value_type>::abs(g)) {
813 c=g/f;
814 r=sqrt((c*c)+1.0);
815 e(i+1)=f*r;
816 c *= (s=1.0/r);
817 } else {
818 s=f/g;
819 r=sqrt((s*s)+1.0);
820 e(i+1)=g*r;
821 s *= (c=1.0/r);
822 }
823 g=d(i+1)-p;
824 r=(d(i)-g)*s+2.0*c*b;
825 p=s*r;
826 d(i+1)=g+p;
827 g=c*r-b;
828 }
829 d(l)=d(l)-p;
830 e(l)=g;
831 e(m)=0.0;
832 }
833 } while (m != l);
834 }
835
836 // order eigenvalues
837 for (i = 0; i < n-1; ++i) {
838 k = i;
839 p = d(i);
840 for (l = i+1; l < n; ++l)
841 if (d(l) < p) {
842 k = l;
843 p = d(l);
844 }
845 d(k) = d(i);
846 d(i) = p;
847 }
848 }
849
850 KOKKOS_INLINE_FUNCTION
851 double
852 Polylib::Serial::
853 GammaFunction(const double x) {
854 double gamma(0);
855
856 if (x == -0.5) gamma = -2.0*sqrt(M_PI);
857 else if (x == 0.0) gamma = 1.0;
858 else if ((x-(ordinal_type)x) == 0.5) {
859 ordinal_type n = (ordinal_type) x;
860 auto tmp = x;
861
862 gamma = sqrt(M_PI);
863 while (n--) {
864 tmp -= 1.0;
865 gamma *= tmp;
866 }
867 } else if ((x-(ordinal_type)x) == 0.0) {
868 ordinal_type n = (ordinal_type) x;
869 auto tmp = x;
870
871 gamma = 1.0;
872 while (--n) {
873 tmp -= 1.0;
874 gamma *= tmp;
875 }
876 } else {
877 INTREPID2_TEST_FOR_ABORT(true,
878 ">>> ERROR (Polylib::Serial): Argument is not of integer or half order.");
879 }
880 return gamma;
881 }
882
883} // end of namespace Intrepid2
884
885#endif
Intrepid2::Util
small utility functions
Definition Intrepid2_Utils.hpp:238
Intrepid2::Polylib::Serial::TriQL
static KOKKOS_INLINE_FUNCTION void TriQL(dViewType d, eViewType e, const ordinal_type n)
QL algorithm for symmetric tridiagonal matrix.
Definition Intrepid2_PolylibDef.hpp:783
Intrepid2::Polylib::Serial::JacobiZeros
static KOKKOS_INLINE_FUNCTION void JacobiZeros(zViewType z, const ordinal_type n, const double alpha, const double beta)
Calculate the n zeros, z, of the Jacobi polynomial, i.e. .
Definition Intrepid2_PolylibDef.hpp:679
Intrepid2::Polylib::Serial::GammaFunction
static KOKKOS_INLINE_FUNCTION double GammaFunction(const double x)
Calculate the Gamma function , , for integer values x and halves.
Definition Intrepid2_PolylibDef.hpp:853
Intrepid2::Polylib::Serial::JacobiPolynomialDerivative
static KOKKOS_INLINE_FUNCTION void JacobiPolynomialDerivative(const ordinal_type np, const zViewType z, polydViewType polyd, const ordinal_type n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition Intrepid2_PolylibDef.hpp:654
Intrepid2::Polylib::Serial::JacobiPolynomial
static KOKKOS_INLINE_FUNCTION void JacobiPolynomial(const ordinal_type np, const zViewType z, polyiViewType poly_in, polydViewType polyd, const ordinal_type n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
Definition Intrepid2_PolylibDef.hpp:573
Generated by doxygen 1.9.8