Intrepid2
Intrepid2_Polylib.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_Polylib.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_HPP__
96#define __INTREPID2_POLYLIB_HPP__
97
98#include "Intrepid2_ConfigDefs.hpp"
99
100#include "Intrepid2_Types.hpp"
101#include "Intrepid2_Utils.hpp"
102
103#include "Kokkos_Core.hpp"
104
105namespace Intrepid2 {
106
204 class Polylib {
205 public:
206
207 static constexpr ordinal_type MaxPolylibIteration = 50;
208 static constexpr ordinal_type MaxPolylibOrder =
211 static constexpr ordinal_type MaxPolylibPoint = MaxPolylibOrder/2+2;
212
213 struct Serial {
214
215 // -----------------------------------------------------------------------
216 // Points and Weights
217 // -----------------------------------------------------------------------
218
236 template<EPolyType polyType>
237 struct Cubature {
238 template<typename zViewType,
239 typename wViewType>
240 KOKKOS_INLINE_FUNCTION
241 static void
242 getValues( zViewType z,
243 wViewType w,
244 const ordinal_type np,
245 const double alpha,
246 const double beta);
247 };
248
249 template<typename zViewType,
250 typename wViewType>
251 KOKKOS_INLINE_FUNCTION
252 static void
253 getCubature( zViewType z,
254 wViewType w,
255 const ordinal_type np,
256 const double alpha,
257 const double beta,
258 const EPolyType poly) {
259 switch (poly) {
260 case POLYTYPE_GAUSS: Polylib::Serial::Cubature<POLYTYPE_GAUSS> ::getValues(z, w, np, alpha, beta); break;
261 case POLYTYPE_GAUSS_RADAU_LEFT: Polylib::Serial::Cubature<POLYTYPE_GAUSS_RADAU_LEFT> ::getValues(z, w, np, alpha, beta); break;
262 case POLYTYPE_GAUSS_RADAU_RIGHT: Polylib::Serial::Cubature<POLYTYPE_GAUSS_RADAU_RIGHT>::getValues(z, w, np, alpha, beta); break;
263 case POLYTYPE_GAUSS_LOBATTO: Polylib::Serial::Cubature<POLYTYPE_GAUSS_LOBATTO> ::getValues(z, w, np, alpha, beta); break;
264 default:
265 INTREPID2_TEST_FOR_ABORT(true,
266 ">>> ERROR (Polylib::Serial::getCubature): Not supported poly type.");
267 break;
268 }
269 }
270
271 // -----------------------------------------------------------------------
272 // Derivative Matrix
273 // -----------------------------------------------------------------------
274
283 template<EPolyType polyType>
284 struct Derivative {
285 template<typename DViewType,
286 typename zViewType>
287 KOKKOS_INLINE_FUNCTION
288 static void
289 getValues( DViewType D,
290 const zViewType z,
291 const ordinal_type np,
292 const double alpha,
293 const double beta);
294 };
295
296 template<typename DViewType,
297 typename zViewType>
298 KOKKOS_INLINE_FUNCTION
299 static void
300 getDerivative( DViewType D,
301 const zViewType z,
302 const ordinal_type np,
303 const double alpha,
304 const double beta,
305 const EPolyType poly) {
306 switch (poly) {
307 case POLYTYPE_GAUSS: Polylib::Serial::Derivative<POLYTYPE_GAUSS> ::getValues(D, z, np, alpha, beta); break;
308 case POLYTYPE_GAUSS_RADAU_LEFT: Polylib::Serial::Derivative<POLYTYPE_GAUSS_RADAU_LEFT> ::getValues(D, z, np, alpha, beta); break;
309 case POLYTYPE_GAUSS_RADAU_RIGHT: Polylib::Serial::Derivative<POLYTYPE_GAUSS_RADAU_RIGHT>::getValues(D, z, np, alpha, beta); break;
310 case POLYTYPE_GAUSS_LOBATTO: Polylib::Serial::Derivative<POLYTYPE_GAUSS_LOBATTO> ::getValues(D, z, np, alpha, beta); break;
311 default:
312 INTREPID2_TEST_FOR_ABORT(true,
313 ">>> ERROR (Polylib::Serial::getDerivative): Not supported poly type.");
314 break;
315 }
316 }
317
318 // -----------------------------------------------------------------------
319 // Lagrangian Interpolants
320 // -----------------------------------------------------------------------
321
381 template<EPolyType polyType>
383 template<typename zViewType>
384 KOKKOS_INLINE_FUNCTION
385 static typename zViewType::value_type
386 getValue(const ordinal_type i,
387 const typename zViewType::value_type z,
388 const zViewType zgj,
389 const ordinal_type np,
390 const double alpha,
391 const double beta);
392 };
393
394 template<typename zViewType>
395 KOKKOS_INLINE_FUNCTION
396 static typename zViewType::value_type
397 getLagrangianInterpolant(const ordinal_type i,
398 const typename zViewType::value_type z,
399 const zViewType zgj,
400 const ordinal_type np,
401 const double alpha,
402 const double beta,
403 const EPolyType poly) {
404 typename zViewType::value_type r_val = 0;
405 switch (poly) {
406 case POLYTYPE_GAUSS: r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS> ::getValue(i, z, zgj, np, alpha, beta); break;
407 case POLYTYPE_GAUSS_RADAU_LEFT: r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_RADAU_LEFT> ::getValue(i, z, zgj, np, alpha, beta); break;
408 case POLYTYPE_GAUSS_RADAU_RIGHT: r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_RADAU_RIGHT>::getValue(i, z, zgj, np, alpha, beta); break;
409 case POLYTYPE_GAUSS_LOBATTO: r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_LOBATTO> ::getValue(i, z, zgj, np, alpha, beta); break;
410 default:
411 INTREPID2_TEST_FOR_ABORT(true,
412 ">>> ERROR (Polylib::Serial::getLagrangianInterpolant): Not supported poly type.");
413 break;
414 }
415 return r_val;
416 }
417
418 // -----------------------------------------------------------------------
419 // Interpolation operators
420 // -----------------------------------------------------------------------
421
432 template<EPolyType polyType>
434 template<typename imViewType,
435 typename zgrjViewType,
436 typename zmViewType>
437 KOKKOS_INLINE_FUNCTION
438 static void
439 getMatrix( imViewType im,
440 const zgrjViewType zgrj,
441 const zmViewType zm,
442 const ordinal_type nz,
443 const ordinal_type mz,
444 const double alpha,
445 const double beta);
446 };
447
448 template<typename imViewType,
449 typename zgrjViewType,
450 typename zmViewType>
451 KOKKOS_INLINE_FUNCTION
452 static void
453 getInterpolationOperator( imViewType im,
454 const zgrjViewType zgrj,
455 const zmViewType zm,
456 const ordinal_type nz,
457 const ordinal_type mz,
458 const double alpha,
459 const double beta,
460 const EPolyType poly) {
461 switch (poly) {
462 case POLYTYPE_GAUSS: Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS> ::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
463 case POLYTYPE_GAUSS_RADAU_LEFT: Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS_RADAU_LEFT> ::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
464 case POLYTYPE_GAUSS_RADAU_RIGHT: Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS_RADAU_RIGHT>::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
465 case POLYTYPE_GAUSS_LOBATTO: Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS_LOBATTO> ::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
466 default:
467 INTREPID2_TEST_FOR_ABORT(true,
468 ">>> ERROR (Polylib::Serial::getInterpolationOperator): Not supported poly type.");
469 break;
470 }
471 }
472
473
474 // -----------------------------------------------------------------------
475 // Polynomial functions
476 // -----------------------------------------------------------------------
477
517 template<typename zViewType,
518 typename polyiViewType,
519 typename polydViewType>
520 KOKKOS_INLINE_FUNCTION
521 static void
522 JacobiPolynomial(const ordinal_type np,
523 const zViewType z,
524 polyiViewType poly_in,
525 polydViewType polyd,
526 const ordinal_type n,
527 const double alpha,
528 const double beta);
529
530
544 template<typename zViewType,
545 typename polydViewType>
546 KOKKOS_INLINE_FUNCTION
547 static void
548 JacobiPolynomialDerivative(const ordinal_type np,
549 const zViewType z,
550 polydViewType polyd,
551 const ordinal_type n,
552 const double alpha,
553 const double beta);
554
555 // -----------------------------------------------------------------------
556 // Helper functions.
557 // -----------------------------------------------------------------------
558
565 template<typename zViewType,
566 bool DeflationEnabled = false>
567 KOKKOS_INLINE_FUNCTION
568 static void
569 JacobiZeros ( zViewType z,
570 const ordinal_type n,
571 const double alpha,
572 const double beta);
573
574 template<typename zViewType>
575 KOKKOS_INLINE_FUNCTION
576 static void
577 JacobiZerosPolyDeflation( zViewType z,
578 const ordinal_type n,
579 const double alpha,
580 const double beta);
581
582 template<typename aViewType>
583 KOKKOS_INLINE_FUNCTION
584 static void
585 JacobiZerosTriDiagonal( aViewType a,
586 const ordinal_type n,
587 const double alpha,
588 const double beta);
589
612 template<typename aViewType,
613 typename bViewType>
614 KOKKOS_INLINE_FUNCTION
615 static void
616 JacobiZeros(aViewType a,
617 bViewType b,
618 const ordinal_type n,
619 const double alpha,
620 const double beta);
621
622
646 template<typename dViewType,
647 typename eViewType>
648 KOKKOS_INLINE_FUNCTION
649 static void
650 TriQL( dViewType d,
651 eViewType e,
652 const ordinal_type n);
653
654
664 KOKKOS_INLINE_FUNCTION
665 static double
666 GammaFunction(const double x);
667
668 };
669
670 // -----------------------------------------------------------------------
671 };
672
673} // end of Intrepid namespace
674
675 // include templated definitions
677
678#endif
Definition file for a set of functions providing orthogonal polynomial calculus and interpolation.
Contains definitions of custom data types in Intrepid2.
Header function for Intrepid2::Util class and other utility functions.
static constexpr ordinal_type MaxCubatureDegreeEdge
The maximum degree of the polynomial that can be integrated exactly by a direct edge rule.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
Providing orthogonal polynomial calculus and interpolation, created by Spencer Sherwin,...
Gauss-Jacobi/Gauss-Radau-Jacobi/Gauss-Lobatto zeros and weights.
Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi/Gauss-Radau-Jacobi/G...
Interpolation Operator from Gauss-Jacobi points to an arbitrary distribution at points zm.
Compute the value of the i th Lagrangian interpolant through the np Gauss-Jacobi/Gauss-Radau-Jacobi/G...
static KOKKOS_INLINE_FUNCTION void TriQL(dViewType d, eViewType e, const ordinal_type n)
QL algorithm for symmetric tridiagonal matrix.
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. .
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...
static KOKKOS_INLINE_FUNCTION double GammaFunction(const double x)
Calculate the Gamma function , , for integer values x and halves.
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.
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, .