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
72namespace 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_PolylibDef.hpp
Definition file for a set of functions providing orthogonal polynomial calculus and interpolation.
Intrepid2_Types.hpp
Contains definitions of custom data types in Intrepid2.
Intrepid2_Utils.hpp
Header function for Intrepid2::Util class and other utility functions.
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::Parameters::MaxOrder
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
Definition Intrepid2_Types.hpp:102
Intrepid2::Polylib
Providing orthogonal polynomial calculus and interpolation, created by Spencer Sherwin,...
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::Derivative
Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi/Gauss-Radau-Jacobi/G...
Definition Intrepid2_Polylib.hpp:252
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::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::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::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::JacobiZeros
static KOKKOS_INLINE_FUNCTION void JacobiZeros(aViewType a, bViewType b, const ordinal_type n, const double alpha, const double beta)
Zero determination through the eigenvalues of a tridiagonal matrix from the three term recursion rela...
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::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
Generated by doxygen 1.9.8