Intrepid2_Polynomials.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//                           Intrepid2 Package
//
// Copyright 2007 NTESS and the Intrepid2 contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef Intrepid2_Polynomials_h
#define Intrepid2_Polynomials_h

#include "Intrepid2_Polylib.hpp"
#include "Intrepid2_Types.hpp"
#include "Intrepid2_Utils.hpp"
#if defined(HAVE_INTREPID2_SACADO) && !defined(KOKKOS_ENABLE_IMPL_VIEW_LEGACY)
#include <Sacado.hpp>
#endif

namespace Intrepid2
{
  namespace Polynomials
  {
    /*
     These polynomials are supplemental to those defined in the Polylib class; there is some overlap.
     We actually take advantage of the overlap in our verification tests, using the Polylib functions
     to verify the ones defined here.  Our interface here is a little simpler, and the functions are a little less
     general than those in Polylib.
     
     We define some polynomial functions that are useful in a variety of contexts.
     In particular, the (integrated) Legendre and Jacobi polynomials are useful in defining
     hierarchical bases.  See in particular:
     
     Federico Fuentes, Brendan Keith, Leszek Demkowicz, Sriram Nagaraj.
     "Orientation embedded high order shape functions for the exact sequence elements of all shapes."
     Computers & Mathematics with Applications, Volume 70, Issue 4, 2015, Pages 353-458, ISSN 0898-1221.
     https://doi.org/10.1016/j.camwa.2015.04.027.
     
     In this implementation, we take care to make minimal assumptions on both the containers
     and the scalar type.  The containers need to support a one-argument operator() for assignment and/or
     lookup (as appropriate).  The scalar type needs to support a cast from Intrepid2::ordinal_type, as well
     as standard arithmetic operations.  In particular, both 1-rank Kokkos::View and Kokkos::DynRankView are
     supported, as are C++ floating point types and Sacado scalar types.
     */
    
    template<typename OutputValueViewType, typename ScalarType>
    KOKKOS_INLINE_FUNCTION void legendreValues(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x)
    {
      if (n >= 0) outputValues(0) = 1.0;
      if (n >= 1) outputValues(1) = x;
      for (int i=2; i<=n; i++)
      {
        const ScalarType i_scalar = ScalarType(i);
        outputValues(i) = (2. - 1. / i_scalar) * x * outputValues(i-1) - (1. - 1. / i_scalar) * outputValues(i-2);
      }
    }
    
    template<typename OutputValueViewType, typename ScalarType>
    KOKKOS_INLINE_FUNCTION void legendreDerivativeValues(OutputValueViewType outputValues, const OutputValueViewType legendreValues, Intrepid2::ordinal_type n, ScalarType x)
    {
      if (n >= 0) outputValues(0) = 0.0;
      if (n >= 1) outputValues(1) = 1.0;
      for (int i=2; i<=n; i++)
      {
        const ScalarType i_scalar = ScalarType(i);
        outputValues(i) = outputValues(i-2) + (2. * i_scalar - 1.) * legendreValues(i-1);
      }
    }
    
    // derivative values can be computed using the Legendre values
    // n: number of Legendre polynomials' derivative values to compute.  outputValues must have at least n+1 entries
    // x: value in [-1, 1]
    // dn: number of derivatives to take.  Must be >= 1.
    template<typename OutputValueViewType, typename PointScalarType>
    KOKKOS_INLINE_FUNCTION void legendreDerivativeValues(OutputValueViewType outputValues, Intrepid2::ordinal_type n, PointScalarType x, Intrepid2::ordinal_type dn)
    {
      const OutputValueViewType nullOutputScalarView;
      const double alpha = 0;
      const double beta = 0;
      const int numPoints = 1;
      
      using Layout = typename NaturalLayoutForType<PointScalarType>::layout;
      
      using UnmanagedPointScalarView = Kokkos::View<PointScalarType*, Layout, Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
      UnmanagedPointScalarView pointView;
#if defined(HAVE_INTREPID2_SACADO) && defined(SACADO_HAS_NEW_KOKKOS_VIEW_IMPL)
      if constexpr (Sacado::is_view_fad<UnmanagedPointScalarView>::value) {
        // This is super weird and doesn't actually wrap the Sacado derivatives
        pointView = UnmanagedPointScalarView(&x.val(), numPoints, 1);
      } else
#endif
      pointView = UnmanagedPointScalarView(&x, numPoints);
      for (int i=0; i<=n; i++)
      {
        auto jacobiValue = Kokkos::subview(outputValues,Kokkos::pair<Intrepid2::ordinal_type,Intrepid2::ordinal_type>(i,i+1));
        jacobiValue(0) = 0.0;
        Intrepid2::Polylib::Serial::JacobiPolynomial(numPoints, pointView, jacobiValue, nullOutputScalarView, i-dn, alpha+dn, beta+dn);
        
        double scaleFactor = 1.0;
        for (int j=1; j<=dn; j++)
        {
          scaleFactor *= 0.5 * (j+alpha+beta+i);
        }
        
        outputValues(i) = jacobiValue(0) * scaleFactor;
      }
    }
    
    template<typename OutputValueViewType, typename ScalarType>
    KOKKOS_INLINE_FUNCTION void shiftedLegendreValues(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x)
    {
      legendreValues(outputValues, n, 2.*x-1.);
    }
    
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledLegendreValues(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      using OutputScalar = typename OutputValueViewType::value_type;
      OutputScalar two_x_minus_t = 2. * x - t;
      OutputScalar t_squared = t * t;
      if (n >= 0) outputValues(0) = 1.0;
      if (n >= 1) outputValues(1) = two_x_minus_t;
      for (int i=2; i<=n; i++)
      {
        const ScalarType one_over_i = 1.0 / ScalarType(i);
        outputValues(i) = one_over_i * ( (2. *i - 1. ) * two_x_minus_t * outputValues(i-1) - (i - 1.) * t_squared * outputValues(i-2));
      }
    }
    
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      // reduced memory version: compute P_i in outputValues
      shiftedScaledLegendreValues(outputValues,n,x,t);
      // keep a copy of the last two P_i values around; update these before overwriting in outputValues
      ScalarType P_i_minus_two, P_i_minus_one;
      if (n >= 0) P_i_minus_two = outputValues(0);
      if (n >= 1) P_i_minus_one = outputValues(1);
      
      if (n >= 0) outputValues(0) = 1.0;
      if (n >= 1) outputValues(1) = x;
      for (int i=2; i<=n; i++)
      {
        const ScalarType & P_i = outputValues(i); // define as P_i just for clarity of the code below
        const ScalarType i_scalar = ScalarType(i);
        ScalarType L_i = (P_i - t * t * P_i_minus_two) /( 2. * (2. * i_scalar - 1.));
        
        // get the next values of P_{i-1} and P_{i-2} before overwriting the P_i value
        P_i_minus_two = P_i_minus_one;
        P_i_minus_one = P_i;
        
        // overwrite P_i value
        outputValues(i) = L_i;
      }
    }
    
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues(OutputValueViewType outputValues, const OutputValueViewType shiftedScaledLegendreValues,
                                                                      Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      // reduced flops version: rely on previously computed P_i
      if (n >= 0) outputValues(0) = 1.0;
      if (n >= 1) outputValues(1) = x;
      for (int i=2; i<=n; i++)
      {
        const ScalarType & P_i           = shiftedScaledLegendreValues(i); // define as P_i just for clarity of the code below
        const ScalarType & P_i_minus_two = shiftedScaledLegendreValues(i-2);
        const ScalarType i_scalar = ScalarType(i);
        outputValues(i) = (P_i - t * t * P_i_minus_two) /( 2. * (2. * i_scalar - 1.));
      }
    }
    
    // the below implementation is commented out for now to guard a certain confusion.
    // the integratedLegendreValues() implementation below is implemented in such a way as to agree, modulo a coordinate transformation, with the
    // shiftedScaledIntegratedLegendreValues() above.  Since the latter really is an integral of shiftedScaledLegendreValues(), the former isn't
    // actually an integral of the (unshifted, unscaled) Legendre polynomials.
//    template<typename OutputValueViewType, typename ScalarType>
//    KOKKOS_INLINE_FUNCTION void integratedLegendreValues(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x)
//    {
//      // to reduce memory requirements, compute P_i in outputValues
//      legendreValues(outputValues,n,x);
//      // keep a copy of the last two P_i values around; update these before overwriting in outputValues
//      ScalarType P_i_minus_two, P_i_minus_one;
//      if (n >= 0) P_i_minus_two = outputValues(0);
//      if (n >= 1) P_i_minus_one = outputValues(1);
//
//      if (n >= 0) outputValues(0) = 1.0;
//      if (n >= 1) outputValues(1) = (x + 1.0) / 2.0;
//      for (int i=2; i<=n; i++)
//      {
//        const ScalarType & P_i = outputValues(i); // define as P_i just for clarity of the code below
//        const ScalarType i_scalar = ScalarType(i);
//        ScalarType L_i = (P_i - P_i_minus_two) /( 2. * (2. * i_scalar - 1.));
//
//        // get the next values of P_{i-1} and P_{i-2} before overwriting the P_i value
//        P_i_minus_two = P_i_minus_one;
//        P_i_minus_one = P_i;
//
//        // overwrite P_i value
//        outputValues(i) = L_i;
//      }
//    }
    
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues_dx(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      if (n >= 0) outputValues(0) = 0.0;
      if (n >= 1) outputValues(1) = 1.0;
      if (n >= 2) outputValues(2) = 2. * x - t;
      for (int i=2; i<=n-1; i++)
      {
        const ScalarType one_over_i = 1.0 / ScalarType(i);
        outputValues(i+1) = one_over_i * (2. * i - 1.) * (2. * x - t) * outputValues(i) -  one_over_i * (i - 1.0) * t * t * outputValues(i-1);
      }
    }
    
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues_dt(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      // memory-conserving version -- place the Legendre values in the final output container
      shiftedScaledLegendreValues(outputValues, n, x, t);
      
      ScalarType P_i_minus_2 = outputValues(0);
      ScalarType P_i_minus_1 = outputValues(1);
      
      if (n >= 0) outputValues(0) = 0.0;
      if (n >= 1) outputValues(1) = 0.0;
      
      for (int i=2; i<=n; i++)
      {
        const ScalarType L_i_dt = -0.5 * (P_i_minus_1 + t * P_i_minus_2);
        
        P_i_minus_2 = P_i_minus_1;
        P_i_minus_1 = outputValues(i);
        
        outputValues(i) = L_i_dt;
      }
    }
    
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues_dt(OutputValueViewType outputValues, const OutputValueViewType shiftedScaledLegendreValues,
                                                                         Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      // reduced flops version: rely on previously computed P_i
      if (n >= 0) outputValues(0) = 0.0;
      if (n >= 1) outputValues(1) = 0.0;
      for (int i=2; i<=n; i++)
      {
        const ScalarType & P_i_minus_1 = shiftedScaledLegendreValues(i-1); // define as P_i just for clarity of the code below
        const ScalarType & P_i_minus_2 = shiftedScaledLegendreValues(i-2);
        outputValues(i) = -0.5 * (P_i_minus_1 + t * P_i_minus_2);
      }
    }
    
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledJacobiValues(OutputValueViewType outputValues, double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      ScalarType two_x_minus_t = 2. * x - t;
      ScalarTypeForScaling alpha_squared_t = alpha * alpha * t;
      
      if (n >= 0) outputValues(0) = 1.0;
      if (n >= 1) outputValues(1) = two_x_minus_t + alpha * x;
      
      for (int i=2; i<=n; i++)
      {
        const ScalarType & P_i_minus_one = outputValues(i-1); // define as P_i just for clarity of the code below
        const ScalarType & P_i_minus_two = outputValues(i-2);
        
        double a_i = (2. * i) * (i + alpha) * (2. * i + alpha - 2.);
        double b_i = 2. * i + alpha - 1.;
        double c_i = (2. * i + alpha) * (2. * i + alpha - 2.);
        double d_i = 2. * (i + alpha - 1.) * (i - 1.) * (2. * i + alpha);
        
        outputValues(i) = (b_i / a_i) * (c_i * two_x_minus_t + alpha_squared_t) * P_i_minus_one - (d_i / a_i) * t * t * P_i_minus_two;
      }
    }
    
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedJacobiValues(OutputValueViewType outputValues, const OutputValueViewType jacobiValues,
                                                                    double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      // reduced flops version: rely on previously computed P_i
      if (n >= 0) outputValues(0) = 1.0;
      if (n >= 1) outputValues(1) = x;
      
      ScalarType t_squared = t * t;
      for (int i=2; i<=n; i++)
      {
        const ScalarType & P_i         = jacobiValues(i);   // define as P_i just for clarity of the code below
        const ScalarType & P_i_minus_1 = jacobiValues(i-1);
        const ScalarType & P_i_minus_2 = jacobiValues(i-2);
        
        double a_i = (i + alpha) / ((2. * i + alpha - 1.) * (2. * i + alpha     ));
        double b_i = alpha       / ((2. * i + alpha - 2.) * (2. * i + alpha     ));
        double c_i = (i - 1.)    / ((2. * i + alpha - 2.) * (2. * i + alpha - 1.));
        
        outputValues(i) = a_i * P_i + b_i * t * P_i_minus_1 - c_i * t_squared * P_i_minus_2;
      }
    }
    
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedJacobiValues(OutputValueViewType outputValues,
                                                                    double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      // memory-conserving version -- place the Jacobi values in the final output container
      shiftedScaledJacobiValues(outputValues, alpha, n, x, t);
      
      ScalarType P_i_minus_2 = outputValues(0);
      ScalarType P_i_minus_1 = outputValues(1);
      
      if (n >= 0) outputValues(0) = 1.0;
      if (n >= 1) outputValues(1) = x;
      
      ScalarType t_squared = t * t;
      for (int i=2; i<=n; i++)
      {
        const ScalarType & P_i = outputValues(i);
        
        double a_i = (i + alpha) / ((2. * i + alpha - 1.) * (2. * i + alpha     ));
        double b_i = alpha       / ((2. * i + alpha - 2.) * (2. * i + alpha     ));
        double c_i = (i - 1.)    / ((2. * i + alpha - 2.) * (2. * i + alpha - 1.));
        
        ScalarType L_i = a_i * P_i + b_i * t * P_i_minus_1 - c_i * t_squared * P_i_minus_2;
        
        P_i_minus_2 = P_i_minus_1;
        P_i_minus_1 = P_i;
        
        outputValues(i) = L_i;
      }
    }
    
    // x derivative of integrated Jacobi is just Jacobi
    // only distinction is in the index -- outputValues indices are shifted by 1 relative to jacobiValues, above
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedJacobiValues_dx(OutputValueViewType outputValues,
                                                                       double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      // rather than repeating the somewhat involved implementation of jacobiValues here,
      // call with (n-1), and then move values accordingly
      shiftedScaledJacobiValues(outputValues, alpha, n-1, x, t);
      
      // forward implementation
      ScalarType nextValue = 0.0;
      ScalarType nextNextValue = 0.0;
      for (int i=0; i<=n-1; i++)
      {
        nextNextValue = outputValues(i);
        outputValues(i) = nextValue;
        nextValue = nextNextValue;
      }
      outputValues(n-1) = nextValue;
    }
    
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedJacobiValues_dt(OutputValueViewType outputValues, const OutputValueViewType jacobiValues,
                                                                       double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      // reduced flops version: rely on previously computed P_i
      if (n >= 0) outputValues(0) = 0.0;
      if (n >= 1) outputValues(1) = 0.0;
      for (int i=2; i<=n; i++)
      {
        const ScalarType & P_i_minus_1 = jacobiValues(i-1); // define as P_i just for clarity of the code below
        const ScalarType & P_i_minus_2 = jacobiValues(i-2);
        outputValues(i) = - (i-1.) / (2. * i - 2. + alpha) * (P_i_minus_1 + t * P_i_minus_2);
      }
    }
    
    template<typename OutputValueViewType, typename ScalarType, typename ScalarTypeForScaling>
    KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedJacobiValues_dt(OutputValueViewType outputValues,
                                                                       double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t)
    {
      // memory-conserving version -- place the Jacobi values in the final output container
      shiftedScaledJacobiValues(outputValues, alpha, n, x, t);
      
      ScalarType P_i_minus_2 = outputValues(0);
      ScalarType P_i_minus_1 = outputValues(1);
      
      if (n >= 0) outputValues(0) = 0.0;
      if (n >= 1) outputValues(1) = 0.0;
      
      for (int i=2; i<=n; i++)
      {
        const ScalarType L_i_dt =  - (i-1.) / (2. * i - 2. + alpha) * (P_i_minus_1 + t * P_i_minus_2);
        
        P_i_minus_2 = P_i_minus_1;
        P_i_minus_1 = outputValues(i);
        
        outputValues(i) = L_i_dt;
      }
    }
  } // namespace Polynomials
} // namespace Intrepid2

#endif /* Intrepid2_Polynomials_h */