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 namespace Intrepid2 {
99 
100  // -----------------------------------------------------------------------
101  // Points and Weights
102  // -----------------------------------------------------------------------
103 
104  template<>
105  template<typename zViewType,
106  typename wViewType>
107  KOKKOS_INLINE_FUNCTION
108  void
109  Polylib::Serial::Cubature<POLYTYPE_GAUSS>::
110  getValues( zViewType z,
111  wViewType w,
112  const ordinal_type np,
113  const double alpha,
114  const double beta) {
115  const double one = 1.0, two = 2.0, apb = alpha + beta;
116  double fac;
117 
118  JacobiZeros(z, np, alpha, beta);
119  JacobiPolynomialDerivative(np, z, w, np, alpha, beta);
120 
121  fac = pow(two, apb + one)*GammaFunction(alpha + np + one)*GammaFunction(beta + np + one);
122  fac /= GammaFunction((np + one))*GammaFunction(apb + np + one);
123 
124  for (ordinal_type i = 0; i < np; ++i)
125  w(i) = fac/(w(i)*w(i)*(one-z(i)*z(i)));
126  }
127 
128  template<>
129  template<typename zViewType,
130  typename wViewType>
131  KOKKOS_INLINE_FUNCTION
132  void
133  Polylib::Serial::Cubature<POLYTYPE_GAUSS_RADAU_LEFT>::
134  getValues( zViewType z,
135  wViewType w,
136  const ordinal_type np,
137  const double alpha,
138  const double beta) {
139  if (np == 1) {
140  z(0) = 0.0;
141  w(0) = 2.0;
142  } else {
143  const double one = 1.0, two = 2.0, apb = alpha + beta;
144  double fac;
145 
146  z(0) = -one;
147 
148  auto z_plus_1 = Kokkos::subview(z, Kokkos::pair<ordinal_type,ordinal_type>(1, z.extent(0)));
149  JacobiZeros(z_plus_1, np-1, alpha, beta+1);
150 
151  Kokkos::View<typename zViewType::value_type*,Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged> null;
152  JacobiPolynomial(np, z, w, null, np-1, alpha, beta);
153 
154  fac = pow(two, apb)*GammaFunction(alpha + np)*GammaFunction(beta + np);
155  fac /= GammaFunction(np)*(beta + np)*GammaFunction(apb + np + 1);
156 
157  for (ordinal_type i = 0; i < np; ++i)
158  w(i) = fac*(1-z(i))/(w(i)*w(i));
159  w(0) *= (beta + one);
160  }
161  }
162 
163  template<>
164  template<typename zViewType,
165  typename wViewType>
166  KOKKOS_INLINE_FUNCTION
167  void
168  Polylib::Serial::Cubature<POLYTYPE_GAUSS_RADAU_RIGHT>::
169  getValues( zViewType z,
170  wViewType w,
171  const ordinal_type np,
172  const double alpha,
173  const double beta) {
174  if (np == 1) {
175  z(0) = 0.0;
176  w(0) = 2.0;
177  } else {
178  const double one = 1.0, two = 2.0, apb = alpha + beta;
179  double fac;
180 
181  JacobiZeros(z, np-1, alpha+1, beta);
182  z(np-1) = one;
183 
184  Kokkos::View<typename zViewType::value_type*,Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged> null;
185  JacobiPolynomial(np, z, w, null, np-1, alpha, beta);
186 
187  fac = pow(two,apb)*GammaFunction(alpha + np)*GammaFunction(beta + np);
188  fac /= GammaFunction(np)*(alpha + np)*GammaFunction(apb + np + 1);
189 
190  for (ordinal_type i = 0; i < np; ++i)
191  w(i) = fac*(1+z(i))/(w(i)*w(i));
192  w(np-1) *= (alpha + one);
193  }
194  }
195 
196  template<>
197  template<typename zViewType,
198  typename wViewType>
199  KOKKOS_INLINE_FUNCTION
200  void
201  Polylib::Serial::Cubature<POLYTYPE_GAUSS_LOBATTO>::
202  getValues( zViewType z,
203  wViewType w,
204  const ordinal_type np,
205  const double alpha,
206  const double beta) {
207  if ( np == 1 ) {
208  z(0) = 0.0;
209  w(0) = 2.0;
210  } else {
211  const double one = 1.0, apb = alpha + beta, two = 2.0;
212  double fac;
213 
214  z(0) = -one;
215  z(np-1) = one;
216 
217  auto z_plus_1 = Kokkos::subview(z, Kokkos::pair<ordinal_type,ordinal_type>(1, z.extent(0)));
218  JacobiZeros(z_plus_1, np-2, alpha+one, beta+one);
219 
220  Kokkos::View<typename zViewType::value_type*,Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged> null;
221  JacobiPolynomial(np, z, w, null, np-1, alpha, beta);
222 
223  fac = pow(two, apb + 1)*GammaFunction(alpha + np)*GammaFunction(beta + np);
224  fac /= (np-1)*GammaFunction(np)*GammaFunction(alpha + beta + np + one);
225 
226  for (ordinal_type i = 0; i < np; ++i)
227  w(i) = fac/(w(i)*w(i));
228 
229  w(0) *= (beta + one);
230  w(np-1) *= (alpha + one);
231  }
232  }
233 
234  // -----------------------------------------------------------------------
235  // Derivatives
236  // -----------------------------------------------------------------------
237 
238  template<>
239  template<typename DViewType,
240  typename zViewType>
241  KOKKOS_INLINE_FUNCTION
242  void
243  Polylib::Serial::Derivative<POLYTYPE_GAUSS>::
244  getValues( DViewType D,
245  const zViewType z,
246  const ordinal_type np,
247  const double alpha,
248  const double beta) {
249  if (np <= 0) {
250  D(0,0) = 0.0;
251  } else {
252  const double one = 1.0, two = 2.0;
253 
254  typename zViewType::value_type pd_buf[MaxPolylibPoint];
255  Kokkos::View<typename zViewType::value_type*,
256  Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged>
257  pd((typename zViewType::pointer_type)&pd_buf[0], MaxPolylibPoint);
258 
259  JacobiPolynomialDerivative(np, z, pd, np, alpha, beta);
260 
261  for (ordinal_type i = 0; i < np; ++i)
262  for (ordinal_type j = 0; j < np; ++j)
263  if (i != j)
264  //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.
265  D(i,j) = pd(i)/(pd(j)*(z(i)-z(j)));
266  else
267  D(i,j) = (alpha - beta + (alpha + beta + two)*z(j))/
268  (two*(one - z(j)*z(j)));
269  }
270  }
271 
272  template<>
273  template<typename DViewType,
274  typename zViewType>
275  KOKKOS_INLINE_FUNCTION
276  void
277  Polylib::Serial::Derivative<POLYTYPE_GAUSS_RADAU_LEFT>::
278  getValues( DViewType D,
279  const zViewType z,
280  const ordinal_type np,
281  const double alpha,
282  const double beta) {
283  if (np <= 0) {
284  D(0,0) = 0.0;
285  } else {
286  const double one = 1.0, two = 2.0;
287 
288  typename zViewType::value_type pd_buf[MaxPolylibPoint];
289  Kokkos::View<typename zViewType::value_type*,
290  Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged>
291  pd((typename zViewType::pointer_type)&pd_buf[0], MaxPolylibPoint);
292 
293  pd(0) = pow(-one,np-1)*GammaFunction(np+beta+one);
294  pd(0) /= GammaFunction(np)*GammaFunction(beta+two);
295 
296  auto pd_plus_1 = Kokkos::subview(pd, Kokkos::pair<ordinal_type,ordinal_type>(1, pd.extent(0)));
297  auto z_plus_1 = Kokkos::subview( z, Kokkos::pair<ordinal_type,ordinal_type>(1, z.extent(0)));
298 
299  JacobiPolynomialDerivative(np-1, z_plus_1, pd_plus_1, np-1, alpha, beta+1);
300  for(ordinal_type i = 1; i < np; ++i)
301  pd(i) *= (1+z(i));
302 
303  for (ordinal_type i = 0; i < np; ++i)
304  for (ordinal_type j = 0; j < np; ++j)
305  if (i != j)
306  D(i,j) = pd(i)/(pd(j)*(z(i)-z(j)));
307  else
308  if (j == 0)
309  D(i,j) = -(np + alpha + beta + one)*(np - one)/
310  (two*(beta + two));
311  else
312  D(i,j) = (alpha - beta + one + (alpha + beta + one)*z(j))/
313  (two*(one - z(j)*z(j)));
314  }
315  }
316 
317  template<>
318  template<typename DViewType,
319  typename zViewType>
320  KOKKOS_INLINE_FUNCTION
321  void
322  Polylib::Serial::Derivative<POLYTYPE_GAUSS_RADAU_RIGHT>::
323  getValues( DViewType D,
324  const zViewType z,
325  const ordinal_type np,
326  const double alpha,
327  const double beta) {
328  if (np <= 0) {
329  D(0,0) = 0.0;
330  } else {
331  const double one = 1.0, two = 2.0;
332 
333  typename zViewType::value_type pd_buf[MaxPolylibPoint];
334  Kokkos::View<typename zViewType::value_type*,
335  Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged>
336  pd((typename zViewType::pointer_type)&pd_buf[0], MaxPolylibPoint);
337 
338  JacobiPolynomialDerivative(np-1, z, pd, np-1, alpha+1, beta);
339  for (ordinal_type i = 0; i < np-1; ++i)
340  pd(i) *= (1-z(i));
341 
342  pd(np-1) = -GammaFunction(np+alpha+one);
343  pd(np-1) /= GammaFunction(np)*GammaFunction(alpha+two);
344 
345  for (ordinal_type i = 0; i < np; ++i)
346  for (ordinal_type j = 0; j < np; ++j)
347  if (i != j)
348  D(i,j) = pd(i)/(pd(j)*(z(i)-z(j)));
349  else
350  if (j == np-1)
351  D(i,j) = (np + alpha + beta + one)*(np - one)/
352  (two*(alpha + two));
353  else
354  D(i,j) = (alpha - beta - one + (alpha + beta + one)*z(j))/
355  (two*(one - z(j)*z(j)));
356  }
357  }
358 
359  template<>
360  template<typename DViewType,
361  typename zViewType>
362  KOKKOS_INLINE_FUNCTION
363  void
364  Polylib::Serial::Derivative<POLYTYPE_GAUSS_LOBATTO>::
365  getValues( DViewType D,
366  const zViewType z,
367  const ordinal_type np,
368  const double alpha,
369  const double beta) {
370  if (np <= 0) {
371  D(0,0) = 0.0;
372  } else {
373  const double one = 1.0, two = 2.0;
374 
375  typename zViewType::value_type pd_buf[MaxPolylibPoint];
376  Kokkos::View<typename zViewType::value_type*,
377  Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged>
378  pd((typename zViewType::pointer_type)&pd_buf[0], MaxPolylibPoint);
379 
380  pd(0) = two*pow(-one,np)*GammaFunction(np + beta);
381  pd(0) /= GammaFunction(np - one)*GammaFunction(beta + two);
382 
383  auto pd_plus_1 = Kokkos::subview(pd, Kokkos::pair<ordinal_type,ordinal_type>(1, pd.extent(0)));
384  auto z_plus_1 = Kokkos::subview( z, Kokkos::pair<ordinal_type,ordinal_type>(1, z.extent(0)));
385 
386  JacobiPolynomialDerivative(np-2, z_plus_1, pd_plus_1, np-2, alpha+1, beta+1);
387  for (ordinal_type i = 1; i < np-1; ++i)
388  pd(i) *= (one-z(i)*z(i));
389 
390  pd(np-1) = -two*GammaFunction(np + alpha);
391  pd(np-1) /= GammaFunction(np - one)*GammaFunction(alpha + two);
392 
393  for (ordinal_type i = 0; i < np; ++i)
394  for (ordinal_type j = 0; j < np; ++j)
395  if (i != j)
396  D(i,j) = pd(i)/(pd(j)*(z(i)-z(j)));
397  else
398  if (j == 0)
399  D(i,j) = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
400  else if (j == np-1)
401  D(i,j) =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
402  else
403  D(i,j) = (alpha - beta + (alpha + beta)*z(j))/
404  (two*(one - z(j)*z(j)));
405  }
406  }
407 
408  // -----------------------------------------------------------------------
409  // Lagrangian Interpolants
410  // -----------------------------------------------------------------------
411 
412  template<>
413  template<typename zViewType>
414  KOKKOS_INLINE_FUNCTION
415  typename zViewType::value_type
416  Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS>::
417  getValue(const ordinal_type i,
418  const typename zViewType::value_type z,
419  const zViewType zgj,
420  const ordinal_type np,
421  const double alpha,
422  const double beta) {
423  const double tol = tolerence();
424 
425  typedef typename zViewType::value_type value_type;
426  typedef typename zViewType::pointer_type pointer_type;
427 
428  value_type h, p_buf, pd_buf, zi_buf = zgj(i);
429  Kokkos::View<value_type*,Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged>
430  p((pointer_type)&p_buf, 1), pd((pointer_type)&pd_buf, 1), zi((pointer_type)&zi_buf, 1),
431  zv(const_cast<pointer_type>(&z), 1), null;
432 
433  const auto dz = z - zi(0);
434  if (Util<value_type>::abs(dz) < tol)
435  return 1.0;
436 
437  JacobiPolynomialDerivative(1, zi, pd, np, alpha, beta);
438  JacobiPolynomial(1, zv, p, null , np, alpha, beta);
439 
440  h = p(0)/(pd(0)*dz);
441 
442  return h;
443  }
444 
445  template<>
446  template<typename zViewType>
447  KOKKOS_INLINE_FUNCTION
448  typename zViewType::value_type
449  Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_RADAU_LEFT>::
450  getValue(const ordinal_type i,
451  const typename zViewType::value_type z,
452  const zViewType zgrj,
453  const ordinal_type np,
454  const double alpha,
455  const double beta) {
456  const double tol = tolerence();
457 
458  typedef typename zViewType::value_type value_type;
459  typedef typename zViewType::pointer_type pointer_type;
460 
461  value_type h, p_buf, pd_buf, zi_buf = zgrj(i);
462  Kokkos::View<value_type*,Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged>
463  p((pointer_type)&p_buf, 1), pd((pointer_type)&pd_buf, 1), zi((pointer_type)&zi_buf, 1),
464  zv(const_cast<pointer_type>(&z), 1), null;
465 
466  const auto dz = z - zi(0);
467  if (Util<value_type>::abs(dz) < tol)
468  return 1.0;
469 
470  JacobiPolynomial(1, zi, p , null, np-1, alpha, beta + 1);
471 
472  // need to use this routine in case zi = -1 or 1
473  JacobiPolynomialDerivative(1, zi, pd, np-1, alpha, beta + 1);
474  h = (1.0 + zi(0))*pd(0) + p(0);
475 
476  JacobiPolynomial(1, zv, p, null, np-1, alpha, beta + 1);
477  h = (1.0 + z)*p(0)/(h*dz);
478 
479  return h;
480  }
481 
482 
483  template<>
484  template<typename zViewType>
485  KOKKOS_INLINE_FUNCTION
486  typename zViewType::value_type
487  Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_RADAU_RIGHT>::
488  getValue(const ordinal_type i,
489  const typename zViewType::value_type z,
490  const zViewType zgrj,
491  const ordinal_type np,
492  const double alpha,
493  const double beta) {
494  const double tol = tolerence();
495 
496  typedef typename zViewType::value_type value_type;
497  typedef typename zViewType::pointer_type pointer_type;
498 
499  value_type h, p_buf, pd_buf, zi_buf = zgrj(i);
500  Kokkos::View<value_type*,Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged>
501  p((pointer_type)&p_buf, 1), pd((pointer_type)&pd_buf, 1), zi((pointer_type)&zi_buf, 1),
502  zv(const_cast<pointer_type>(&z), 1), null;
503 
504  const auto dz = z - zi(0);
505  if (Util<value_type>::abs(dz) < tol)
506  return 1.0;
507 
508  JacobiPolynomial(1, zi, p , null, np-1, alpha+1, beta);
509 
510  // need to use this routine in case z = -1 or 1
511  JacobiPolynomialDerivative(1, zi, pd, np-1, alpha+1, beta);
512  h = (1.0 - zi(0))*pd(0) - p(0);
513 
514  JacobiPolynomial (1, zv, p, null, np-1, alpha+1, beta);
515  h = (1.0 - z)*p(0)/(h*dz);
516 
517  return h;
518  }
519 
520 
521  template<>
522  template<typename zViewType>
523  KOKKOS_INLINE_FUNCTION
524  typename zViewType::value_type
525  Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_LOBATTO>::
526  getValue(const ordinal_type i,
527  const typename zViewType::value_type z,
528  const zViewType zglj,
529  const ordinal_type np,
530  const double alpha,
531  const double beta) {
532  const double one = 1.0, two = 2.0, tol = tolerence();
533 
534  typedef typename zViewType::value_type value_type;
535  typedef typename zViewType::pointer_type pointer_type;
536 
537  value_type h, p_buf, pd_buf, zi_buf = zglj(i);
538  Kokkos::View<value_type*,Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged>
539  p((pointer_type)&p_buf, 1), pd((pointer_type)&pd_buf, 1), zi((pointer_type)&zi_buf, 1),
540  zv(const_cast<pointer_type>(&z), 1), null;
541 
542  const auto dz = z - zi(0);
543  if (Util<value_type>::abs(dz) < tol)
544  return 1.0;
545 
546  JacobiPolynomial(1, zi, p , null, np-2, alpha + one, beta + one);
547 
548  // need to use this routine in case z = -1 or 1
549  JacobiPolynomialDerivative(1, zi, pd, np-2, alpha + one, beta + one);
550  h = (one - zi(0)*zi(0))*pd(0) - two*zi(0)*p(0);
551 
552  JacobiPolynomial(1, zv, p, null, np-2, alpha + one, beta + one);
553  h = (one - z*z)*p(0)/(h*dz);
554 
555  return h;
556  }
557 
558  // -----------------------------------------------------------------------
559  // Interpolation Operator
560  // -----------------------------------------------------------------------
561 
562  template<EPolyType polyType>
563  template<typename imViewType,
564  typename zgrjViewType,
565  typename zmViewType>
566  KOKKOS_INLINE_FUNCTION
567  void
568  Polylib::Serial::InterpolationOperator<polyType>::
569  getMatrix( imViewType im,
570  const zgrjViewType zgrj,
571  const zmViewType zm,
572  const ordinal_type nz,
573  const ordinal_type mz,
574  const double alpha,
575  const double beta) {
576  for (ordinal_type i = 0; i < mz; ++i) {
577  const auto zp = zm(i);
578  for (ordinal_type j = 0; j < nz; ++j)
579  im(i, j) = LagrangianInterpolant<polyType>::getValue(j, zp, zgrj, nz, alpha, beta);
580  }
581  }
582 
583  // -----------------------------------------------------------------------
584 
585  template<typename zViewType,
586  typename polyiViewType,
587  typename polydViewType>
588  KOKKOS_INLINE_FUNCTION
589  void
590  Polylib::Serial::
591  JacobiPolynomial(const ordinal_type np,
592  const zViewType z,
593  polyiViewType polyi,
594  polydViewType polyd,
595  const ordinal_type n,
596  const double alpha,
597  const double beta) {
598  const double zero = 0.0, one = 1.0, two = 2.0;
599 
600  if (! np) {
601  return;
602  }
603 
604  if (n == 0) {
605  if (polyi.data())
606  for (ordinal_type i = 0; i < np; ++i)
607  polyi(i) = one;
608  if (polyd.data())
609  for (ordinal_type i = 0; i < np; ++i)
610  polyd(i) = zero;
611  } else if (n == 1) {
612  if (polyi.data())
613  for (ordinal_type i = 0; i < np; ++i)
614  polyi(i) = 0.5*(alpha - beta + (alpha + beta + two)*z(i));
615  if (polyd.data())
616  for (ordinal_type i = 0; i < np; ++i)
617  polyd(i) = 0.5*(alpha + beta + two);
618  } else {
619  double a1, a2, a3, a4;
620  const double apb = alpha + beta;
621 
622  typename polyiViewType::value_type
623  poly[MaxPolylibPoint]={}, polyn1[MaxPolylibPoint]={}, polyn2[MaxPolylibPoint]={};
624 
625  if (polyi.data())
626  for (ordinal_type i=0;i<np;++i)
627  poly[i] = polyi(i);
628 
629  for (ordinal_type i = 0; i < np; ++i) {
630  polyn2[i] = one;
631  polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z(i));
632  }
633 
634  for (auto k = 2; k <= n; ++k) {
635  a1 = two*k*(k + apb)*(two*k + apb - two);
636  a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
637  a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
638  a4 = two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
639 
640  a2 /= a1;
641  a3 /= a1;
642  a4 /= a1;
643 
644  for (ordinal_type i = 0; i < np; ++i) {
645  poly [i] = (a2 + a3*z(i))*polyn1[i] - a4*polyn2[i];
646  polyn2[i] = polyn1[i];
647  polyn1[i] = poly [i];
648  }
649  }
650 
651  if (polyd.data()) {
652  a1 = n*(alpha - beta);
653  a2 = n*(two*n + alpha + beta);
654  a3 = two*(n + alpha)*(n + beta);
655  a4 = (two*n + alpha + beta);
656  a1 /= a4;
657  a2 /= a4;
658  a3 /= a4;
659 
660  // note polyn2 points to polyn1 at end of poly iterations
661  for (ordinal_type i = 0; i < np; ++i) {
662  polyd(i) = (a1- a2*z(i))*poly[i] + a3*polyn2[i];
663  polyd(i) /= (one - z(i)*z(i));
664  }
665  }
666 
667  if (polyi.data())
668  for (ordinal_type i=0;i<np;++i)
669  polyi(i) = poly[i];
670  }
671  }
672 
673  template<typename zViewType,
674  typename polydViewType>
675  KOKKOS_INLINE_FUNCTION
676  void
677  Polylib::Serial::
678  JacobiPolynomialDerivative(const ordinal_type np,
679  const zViewType z,
680  polydViewType polyd,
681  const ordinal_type n,
682  const double alpha,
683  const double beta) {
684  const double one = 1.0;
685  if (n == 0)
686  for(ordinal_type i = 0; i < np; ++i)
687  polyd(i) = 0.0;
688  else {
689  Kokkos::View<typename polydViewType::value_type*,Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged> null;
690  JacobiPolynomial(np, z, polyd, null, n-1, alpha+one, beta+one);
691  for(ordinal_type i = 0; i < np; ++i)
692  polyd(i) *= 0.5*(alpha + beta + n + one);
693  }
694  }
695 
696  // -----------------------------------------------------------------------
697 
698  template<typename zViewType,
699  bool DeflationEnabled>
700  KOKKOS_INLINE_FUNCTION
701  void
702  Polylib::Serial::
703  JacobiZeros( zViewType z,
704  const ordinal_type n,
705  const double alpha,
706  const double beta) {
707  if (DeflationEnabled)
708  JacobiZerosPolyDeflation(z, n, alpha, beta);
709  else
710  JacobiZerosTriDiagonal(z, n, alpha, beta);
711  }
712 
713  template<typename zViewType>
714  KOKKOS_INLINE_FUNCTION
715  void
716  Polylib::Serial::
717  JacobiZerosPolyDeflation( zViewType z,
718  const ordinal_type n,
719  const double alpha,
720  const double beta) {
721  if(!n)
722  return;
723 
724  const double dth = M_PI/(2.0*n);
725  const double one = 1.0, two = 2.0;
726  const double tol = tolerence();
727 
728  typedef typename zViewType::value_type value_type;
729  typedef typename zViewType::pointer_type pointer_type;
730 
731  value_type r_buf, poly_buf, pder_buf;
732  Kokkos::View<value_type*,Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged>
733  poly((pointer_type)&poly_buf, 1), pder((pointer_type)&pder_buf, 1), r((pointer_type)&r_buf, 1);
734 
735  value_type rlast = 0.0;
736  for (auto k = 0; k < n; ++k) {
737  r(0) = -cos((two*(double)k + one) * dth);
738  if (k)
739  r(0) = 0.5*(r(0) + rlast);
740 
741  for (ordinal_type j = 1; j < MaxPolylibIteration; ++j) {
742  JacobiPolynomial(1, r, poly, pder, n, alpha, beta);
743 
744  value_type sum = 0;
745  for (ordinal_type i = 0; i < k; ++i)
746  sum += one/(r(0) - z(i));
747 
748  const value_type delr = -poly(0) / (pder(0) - sum * poly(0));
749  r(0) += delr;
750 
751  if( Util<value_type>::abs(delr) < tol )
752  break;
753  }
754  z(k) = r(0);
755  rlast = r(0);
756  }
757  }
758 
759  template<typename aViewType>
760  KOKKOS_INLINE_FUNCTION
761  void
762  Polylib::Serial::
763  JacobiZerosTriDiagonal( aViewType a,
764  const ordinal_type n,
765  const double alpha,
766  const double beta) {
767  if(!n)
768  return;
769 
770  typedef typename aViewType::value_type value_type;
771  typedef typename aViewType::pointer_type pointer_type;
772 
773  value_type b_buf[MaxPolylibPoint];
774  Kokkos::View<value_type*,Kokkos::Impl::ActiveExecutionMemorySpace,Kokkos::MemoryUnmanaged>
775  b((pointer_type)&b_buf[0], MaxPolylibPoint);
776 
777  // generate normalised terms
778  const double apb = alpha + beta;
779  double apbi = 2.0 + apb;
780 
781  b(n-1) = pow(2.0,apb+1.0)*GammaFunction(alpha+1.0)*GammaFunction(beta+1.0)/GammaFunction(apbi);
782  a(0) = (beta-alpha)/apbi;
783  b(0) = sqrt(4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
784 
785  const double a2b2 = beta*beta-alpha*alpha;
786  for (ordinal_type i = 1; i < n-1; ++i) {
787  apbi = 2.0*(i+1) + apb;
788  a(i) = a2b2/((apbi-2.0)*apbi);
789  b(i) = sqrt(4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
790  ((apbi*apbi-1)*apbi*apbi));
791  }
792 
793  apbi = 2.0*n + apb;
794  //a(n-1) = a2b2/((apbi-2.0)*apbi); // THIS IS A BUG!!!
795  if (n>1) a(n-1) = a2b2/((apbi-2.0)*apbi);
796 
797  // find eigenvalues
798  TriQL(a, b, n);
799  }
800 
801 
802  template<typename dViewType,
803  typename eViewType>
804  KOKKOS_INLINE_FUNCTION
805  void
806  Polylib::Serial::
807  TriQL( dViewType d,
808  eViewType e,
809  const ordinal_type n) {
810  ordinal_type m,l,iter,i,k;
811 
812  typedef typename dViewType::value_type value_type;
813  value_type s,r,p,g,f,dd,c,b;
814 
815  for (l=0; l<n; ++l) {
816  iter=0;
817  do {
818  for (m=l; m<n-1; ++m) {
819  dd=Util<value_type>::abs(d(m))+Util<value_type>::abs(d(m+1));
820  if (Util<value_type>::abs(e(m))+dd == dd) break;
821  }
822  if (m != l) {
823  if (iter++ == MaxPolylibIteration) {
824  INTREPID2_TEST_FOR_ABORT(true,
825  ">>> ERROR (Polylib::Serial): Too many iterations in TQLI.");
826  }
827  g=(d(l+1)-d(l))/(2.0*e(l));
828  r=sqrt((g*g)+1.0);
829  //g=d(m)-d(l)+e(l)/(g+sign(r,g));
830  g=d(m)-d(l)+e(l)/(g+((g)<0 ? value_type(-Util<value_type>::abs(r)) : Util<value_type>::abs(r)));
831  s=c=1.0;
832  p=0.0;
833  for (i=m-1; i>=l; i--) {
834  f=s*e(i);
835  b=c*e(i);
836  if (Util<value_type>::abs(f) >= Util<value_type>::abs(g)) {
837  c=g/f;
838  r=sqrt((c*c)+1.0);
839  e(i+1)=f*r;
840  c *= (s=1.0/r);
841  } else {
842  s=f/g;
843  r=sqrt((s*s)+1.0);
844  e(i+1)=g*r;
845  s *= (c=1.0/r);
846  }
847  g=d(i+1)-p;
848  r=(d(i)-g)*s+2.0*c*b;
849  p=s*r;
850  d(i+1)=g+p;
851  g=c*r-b;
852  }
853  d(l)=d(l)-p;
854  e(l)=g;
855  e(m)=0.0;
856  }
857  } while (m != l);
858  }
859 
860  // order eigenvalues
861  for (i = 0; i < n-1; ++i) {
862  k = i;
863  p = d(i);
864  for (l = i+1; l < n; ++l)
865  if (d(l) < p) {
866  k = l;
867  p = d(l);
868  }
869  d(k) = d(i);
870  d(i) = p;
871  }
872  }
873 
874  KOKKOS_INLINE_FUNCTION
875  double
876  Polylib::Serial::
877  GammaFunction(const double x) {
878  double gamma(0);
879 
880  if (x == -0.5) gamma = -2.0*sqrt(M_PI);
881  else if (x == 0.0) gamma = 1.0;
882  else if ((x-(ordinal_type)x) == 0.5) {
883  ordinal_type n = (ordinal_type) x;
884  auto tmp = x;
885 
886  gamma = sqrt(M_PI);
887  while (n--) {
888  tmp -= 1.0;
889  gamma *= tmp;
890  }
891  } else if ((x-(ordinal_type)x) == 0.0) {
892  ordinal_type n = (ordinal_type) x;
893  auto tmp = x;
894 
895  gamma = 1.0;
896  while (--n) {
897  tmp -= 1.0;
898  gamma *= tmp;
899  }
900  } else {
901  INTREPID2_TEST_FOR_ABORT(true,
902  ">>> ERROR (Polylib::Serial): Argument is not of integer or half order.");
903  }
904  return gamma;
905  }
906 
907 } // end of namespace Intrepid2
908 
909 #endif
Intrepid2::Util
small utility functions
Definition: Intrepid2_Utils.hpp:248
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:678
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:591
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:703
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:877
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:807
Generated by   doxygen 1.8.5