Intrepid2
Intrepid2_Polylib.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_Polylib.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_HPP__
63 #define __INTREPID2_POLYLIB_HPP__
64 
65 #include "Intrepid2_ConfigDefs.hpp"
66 
67 #include "Intrepid2_Types.hpp"
68 #include "Intrepid2_Utils.hpp"
69 
70 #include "Kokkos_Core.hpp"
71 
72 namespace Intrepid2 {
73 
171  class Polylib {
172  public:
173 
174  static constexpr ordinal_type MaxPolylibIteration = 50;
175  static constexpr ordinal_type MaxPolylibOrder =
176  (Parameters::MaxOrder > Parameters::MaxCubatureDegreeEdge) ? Parameters::MaxOrder :
177  Parameters::MaxCubatureDegreeEdge;
178  // NVR: HVOL bases on tri/tet use Polylib with order + spaceDim + 2 points; in 3D this can be up to Parameters::MaxOrder + 5.
179  static constexpr ordinal_type MaxPolylibPoint = (MaxPolylibOrder/2+2 > Parameters::MaxOrder + 5) ? MaxPolylibOrder/2+2 : Parameters::MaxOrder + 5;
180 
181  struct Serial {
182 
183  // -----------------------------------------------------------------------
184  // Points and Weights
185  // -----------------------------------------------------------------------
186 
204  template<EPolyType polyType>
205  struct Cubature {
206  template<typename zViewType,
207  typename wViewType>
208  KOKKOS_INLINE_FUNCTION
209  static void
210  getValues( zViewType z,
211  wViewType w,
212  const ordinal_type np,
213  const double alpha,
214  const double beta);
215  };
216 
217  template<typename zViewType,
218  typename wViewType>
219  KOKKOS_INLINE_FUNCTION
220  static void
221  getCubature( zViewType z,
222  wViewType w,
223  const ordinal_type np,
224  const double alpha,
225  const double beta,
226  const EPolyType poly) {
227  switch (poly) {
228  case POLYTYPE_GAUSS: Polylib::Serial::Cubature<POLYTYPE_GAUSS> ::getValues(z, w, np, alpha, beta); break;
229  case POLYTYPE_GAUSS_RADAU_LEFT: Polylib::Serial::Cubature<POLYTYPE_GAUSS_RADAU_LEFT> ::getValues(z, w, np, alpha, beta); break;
230  case POLYTYPE_GAUSS_RADAU_RIGHT: Polylib::Serial::Cubature<POLYTYPE_GAUSS_RADAU_RIGHT>::getValues(z, w, np, alpha, beta); break;
231  case POLYTYPE_GAUSS_LOBATTO: Polylib::Serial::Cubature<POLYTYPE_GAUSS_LOBATTO> ::getValues(z, w, np, alpha, beta); break;
232  default:
233  INTREPID2_TEST_FOR_ABORT(true,
234  ">>> ERROR (Polylib::Serial::getCubature): Not supported poly type.");
235  break;
236  }
237  }
238 
239  // -----------------------------------------------------------------------
240  // Derivative Matrix
241  // -----------------------------------------------------------------------
242 
251  template<EPolyType polyType>
252  struct Derivative {
253  template<typename DViewType,
254  typename zViewType>
255  KOKKOS_INLINE_FUNCTION
256  static void
257  getValues( DViewType D,
258  const zViewType z,
259  const ordinal_type np,
260  const double alpha,
261  const double beta);
262  };
263 
264  template<typename DViewType,
265  typename zViewType>
266  KOKKOS_INLINE_FUNCTION
267  static void
268  getDerivative( DViewType D,
269  const zViewType z,
270  const ordinal_type np,
271  const double alpha,
272  const double beta,
273  const EPolyType poly) {
274  switch (poly) {
275  case POLYTYPE_GAUSS: Polylib::Serial::Derivative<POLYTYPE_GAUSS> ::getValues(D, z, np, alpha, beta); break;
276  case POLYTYPE_GAUSS_RADAU_LEFT: Polylib::Serial::Derivative<POLYTYPE_GAUSS_RADAU_LEFT> ::getValues(D, z, np, alpha, beta); break;
277  case POLYTYPE_GAUSS_RADAU_RIGHT: Polylib::Serial::Derivative<POLYTYPE_GAUSS_RADAU_RIGHT>::getValues(D, z, np, alpha, beta); break;
278  case POLYTYPE_GAUSS_LOBATTO: Polylib::Serial::Derivative<POLYTYPE_GAUSS_LOBATTO> ::getValues(D, z, np, alpha, beta); break;
279  default:
280  INTREPID2_TEST_FOR_ABORT(true,
281  ">>> ERROR (Polylib::Serial::getDerivative): Not supported poly type.");
282  break;
283  }
284  }
285 
286  // -----------------------------------------------------------------------
287  // Lagrangian Interpolants
288  // -----------------------------------------------------------------------
289 
349  template<EPolyType polyType>
350  struct LagrangianInterpolant {
351  template<typename zViewType>
352  KOKKOS_INLINE_FUNCTION
353  static typename zViewType::value_type
354  getValue(const ordinal_type i,
355  const typename zViewType::value_type z,
356  const zViewType zgj,
357  const ordinal_type np,
358  const double alpha,
359  const double beta);
360  };
361 
362  template<typename zViewType>
363  KOKKOS_INLINE_FUNCTION
364  static typename zViewType::value_type
365  getLagrangianInterpolant(const ordinal_type i,
366  const typename zViewType::value_type z,
367  const zViewType zgj,
368  const ordinal_type np,
369  const double alpha,
370  const double beta,
371  const EPolyType poly) {
372  typename zViewType::value_type r_val = 0;
373  switch (poly) {
374  case POLYTYPE_GAUSS: r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS> ::getValue(i, z, zgj, np, alpha, beta); break;
375  case POLYTYPE_GAUSS_RADAU_LEFT: r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_RADAU_LEFT> ::getValue(i, z, zgj, np, alpha, beta); break;
376  case POLYTYPE_GAUSS_RADAU_RIGHT: r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_RADAU_RIGHT>::getValue(i, z, zgj, np, alpha, beta); break;
377  case POLYTYPE_GAUSS_LOBATTO: r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_LOBATTO> ::getValue(i, z, zgj, np, alpha, beta); break;
378  default:
379  INTREPID2_TEST_FOR_ABORT(true,
380  ">>> ERROR (Polylib::Serial::getLagrangianInterpolant): Not supported poly type.");
381  break;
382  }
383  return r_val;
384  }
385 
386  // -----------------------------------------------------------------------
387  // Interpolation operators
388  // -----------------------------------------------------------------------
389 
400  template<EPolyType polyType>
401  struct InterpolationOperator {
402  template<typename imViewType,
403  typename zgrjViewType,
404  typename zmViewType>
405  KOKKOS_INLINE_FUNCTION
406  static void
407  getMatrix( imViewType im,
408  const zgrjViewType zgrj,
409  const zmViewType zm,
410  const ordinal_type nz,
411  const ordinal_type mz,
412  const double alpha,
413  const double beta);
414  };
415 
416  template<typename imViewType,
417  typename zgrjViewType,
418  typename zmViewType>
419  KOKKOS_INLINE_FUNCTION
420  static void
421  getInterpolationOperator( imViewType im,
422  const zgrjViewType zgrj,
423  const zmViewType zm,
424  const ordinal_type nz,
425  const ordinal_type mz,
426  const double alpha,
427  const double beta,
428  const EPolyType poly) {
429  switch (poly) {
430  case POLYTYPE_GAUSS: Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS> ::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
431  case POLYTYPE_GAUSS_RADAU_LEFT: Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS_RADAU_LEFT> ::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
432  case POLYTYPE_GAUSS_RADAU_RIGHT: Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS_RADAU_RIGHT>::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
433  case POLYTYPE_GAUSS_LOBATTO: Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS_LOBATTO> ::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
434  default:
435  INTREPID2_TEST_FOR_ABORT(true,
436  ">>> ERROR (Polylib::Serial::getInterpolationOperator): Not supported poly type.");
437  break;
438  }
439  }
440 
441 
442  // -----------------------------------------------------------------------
443  // Polynomial functions
444  // -----------------------------------------------------------------------
445 
485  template<typename zViewType,
486  typename polyiViewType,
487  typename polydViewType>
488  KOKKOS_INLINE_FUNCTION
489  static void
490  JacobiPolynomial(const ordinal_type np,
491  const zViewType z,
492  polyiViewType poly_in,
493  polydViewType polyd,
494  const ordinal_type n,
495  const double alpha,
496  const double beta);
497 
498 
512  template<typename zViewType,
513  typename polydViewType>
514  KOKKOS_INLINE_FUNCTION
515  static void
516  JacobiPolynomialDerivative(const ordinal_type np,
517  const zViewType z,
518  polydViewType polyd,
519  const ordinal_type n,
520  const double alpha,
521  const double beta);
522 
523  // -----------------------------------------------------------------------
524  // Helper functions.
525  // -----------------------------------------------------------------------
526 
533  template<typename zViewType,
534  bool DeflationEnabled = false>
535  KOKKOS_INLINE_FUNCTION
536  static void
537  JacobiZeros ( zViewType z,
538  const ordinal_type n,
539  const double alpha,
540  const double beta);
541 
542  template<typename zViewType>
543  KOKKOS_INLINE_FUNCTION
544  static void
545  JacobiZerosPolyDeflation( zViewType z,
546  const ordinal_type n,
547  const double alpha,
548  const double beta);
549 
550  template<typename aViewType>
551  KOKKOS_INLINE_FUNCTION
552  static void
553  JacobiZerosTriDiagonal( aViewType a,
554  const ordinal_type n,
555  const double alpha,
556  const double beta);
557 
580  template<typename aViewType,
581  typename bViewType>
582  KOKKOS_INLINE_FUNCTION
583  static void
584  JacobiZeros(aViewType a,
585  bViewType b,
586  const ordinal_type n,
587  const double alpha,
588  const double beta);
589 
590 
614  template<typename dViewType,
615  typename eViewType>
616  KOKKOS_INLINE_FUNCTION
617  static void
618  TriQL( dViewType d,
619  eViewType e,
620  const ordinal_type n);
621 
622 
632  KOKKOS_INLINE_FUNCTION
633  static double
634  GammaFunction(const double x);
635 
636  };
637 
638  // -----------------------------------------------------------------------
639  };
640 
641 } // end of Intrepid namespace
642 
643  // include templated definitions
644 #include <Intrepid2_PolylibDef.hpp>
645 
646 #endif
Intrepid2::Polylib
Providing orthogonal polynomial calculus and interpolation, created by Spencer Sherwin, Aeronautics, Imperial College London, modified and redistributed by D. Ridzal.
Definition: Intrepid2_Polylib.hpp:171
Intrepid2::Polylib::Serial::Cubature
Gauss-Jacobi/Gauss-Radau-Jacobi/Gauss-Lobatto zeros and weights.
Definition: Intrepid2_Polylib.hpp:205
Intrepid2::Polylib::Serial::LagrangianInterpolant
Compute the value of the i th Lagrangian interpolant through the np Gauss-Jacobi/Gauss-Radau-Jacobi/G...
Definition: Intrepid2_Polylib.hpp:350
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_Utils.hpp
Header function for Intrepid2::Util class and other utility functions.
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_PolylibDef.hpp
Definition file for a set of functions providing orthogonal polynomial calculus and interpolation...
Intrepid2::Polylib::Serial::InterpolationOperator
Interpolation Operator from Gauss-Jacobi points to an arbitrary distribution at points zm...
Definition: Intrepid2_Polylib.hpp:401
Intrepid2::Polylib::Serial::Derivative
Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi/Gauss-Radau-Jacobi/G...
Definition: Intrepid2_Polylib.hpp:252
Intrepid2_Types.hpp
Contains definitions of custom data types in Intrepid2.
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::Parameters::MaxCubatureDegreeEdge
static constexpr ordinal_type MaxCubatureDegreeEdge
The maximum degree of the polynomial that can be integrated exactly by a direct edge rule...
Definition: Intrepid2_Types.hpp:106
Intrepid2::Polylib::Serial::TriQL
static KOKKOS_INLINE_FUNCTION void TriQL(dViewType d, eViewType e, const ordinal_type n)
QL algorithm for symmetric tridiagonal matrix.
Definition: Intrepid2_PolylibDef.hpp:783
Intrepid2::Parameters::MaxOrder
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
Definition: Intrepid2_Types.hpp:102
Generated by   doxygen 1.8.5