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 
73 namespace 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:237
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
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::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
Generated by   doxygen 1.8.5