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 
106 namespace 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
598  Polylib::Serial::
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
685  Polylib::Serial::
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
710  Polylib::Serial::
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
814  Polylib::Serial::
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) {
827  dd=Util<value_type>::abs(d(m))+Util<value_type>::abs(d(m+1));
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);
844  if (Util<value_type>::abs(f) >= Util<value_type>::abs(g)) {
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
884  Polylib::Serial::
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
Intrepid2::Util
small utility functions
Definition: Intrepid2_Utils.hpp:270
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:686
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:599
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:711
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:885
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:815
Generated by   doxygen 1.8.5