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