Intrepid2_Polylib.hpp

Go to the documentation of this file.

/*
// @HEADER
// ************************************************************************
//
//                           Intrepid2 Package
//                 Copyright (2007) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Kyungjoo Kim  (kyukim@sandia.gov), or
//                    Mauro Perego  (mperego@sandia.gov)
//
// ************************************************************************
// @HEADER
*/

//
// File: Intrepid_Polylib.hpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description:
// This file is redistributed with the Intrepid package. It should be used
// in accordance with the above MIT license, at the request of the authors.
// This file is NOT covered by the usual Intrepid/Trilinos LGPL license.
//
// Origin: Nektar++ library, http://www.nektar.info, downloaded on
//         March 10, 2009.
//


#ifndef __INTREPID2_POLYLIB_HPP__
#define __INTREPID2_POLYLIB_HPP__

#include "Intrepid2_ConfigDefs.hpp"

#include "Intrepid2_Types.hpp"
#include "Intrepid2_Utils.hpp"

#include "Kokkos_Core.hpp"

namespace Intrepid2 {

  class Polylib {
  public:

    static constexpr ordinal_type  MaxPolylibIteration = 50;
    static constexpr ordinal_type  MaxPolylibOrder =
        (Parameters::MaxOrder > Parameters::MaxCubatureDegreeEdge) ? Parameters::MaxOrder :
                                                                     Parameters::MaxCubatureDegreeEdge;
    static constexpr ordinal_type  MaxPolylibPoint = MaxPolylibOrder/2+2;

    struct Serial {

      // -----------------------------------------------------------------------
      // Points and Weights
      // -----------------------------------------------------------------------
      
      template<EPolyType polyType>
      struct Cubature {
        template<typename zViewType,
                 typename wViewType>
        KOKKOS_INLINE_FUNCTION
        static void
        getValues(      zViewType z, 
                        wViewType w, 
                  const ordinal_type np, 
                  const double alpha, 
                  const double beta);
      };

      template<typename zViewType,
               typename wViewType>
      KOKKOS_INLINE_FUNCTION
      static void
      getCubature(      zViewType z, 
                        wViewType w, 
                  const ordinal_type np, 
                  const double alpha, 
                  const double beta,
                  const EPolyType poly) {
        switch (poly) {
        case POLYTYPE_GAUSS:             Polylib::Serial::Cubature<POLYTYPE_GAUSS>            ::getValues(z, w, np, alpha, beta); break;
        case POLYTYPE_GAUSS_RADAU_LEFT:  Polylib::Serial::Cubature<POLYTYPE_GAUSS_RADAU_LEFT> ::getValues(z, w, np, alpha, beta); break;
        case POLYTYPE_GAUSS_RADAU_RIGHT: Polylib::Serial::Cubature<POLYTYPE_GAUSS_RADAU_RIGHT>::getValues(z, w, np, alpha, beta); break;
        case POLYTYPE_GAUSS_LOBATTO:     Polylib::Serial::Cubature<POLYTYPE_GAUSS_LOBATTO>    ::getValues(z, w, np, alpha, beta); break;
        default:
          INTREPID2_TEST_FOR_ABORT(true,
                                   ">>> ERROR (Polylib::Serial::getCubature): Not supported poly type.");
          break;
        }
      }
     
      // -----------------------------------------------------------------------
      // Derivative Matrix
      // -----------------------------------------------------------------------

      template<EPolyType polyType>
      struct Derivative {
        template<typename DViewType,
                 typename zViewType>
        KOKKOS_INLINE_FUNCTION
        static void
        getValues(      DViewType D,  
                  const zViewType z, 
                  const ordinal_type np, 
                  const double alpha, 
                  const double beta);
      };

      template<typename DViewType,
               typename zViewType>
      KOKKOS_INLINE_FUNCTION
      static void
      getDerivative(      DViewType D,  
                    const zViewType z, 
                    const ordinal_type np, 
                    const double alpha, 
                    const double beta,
                    const EPolyType poly) {
        switch (poly) {
        case POLYTYPE_GAUSS:             Polylib::Serial::Derivative<POLYTYPE_GAUSS>            ::getValues(D, z, np, alpha, beta); break;
        case POLYTYPE_GAUSS_RADAU_LEFT:  Polylib::Serial::Derivative<POLYTYPE_GAUSS_RADAU_LEFT> ::getValues(D, z, np, alpha, beta); break;
        case POLYTYPE_GAUSS_RADAU_RIGHT: Polylib::Serial::Derivative<POLYTYPE_GAUSS_RADAU_RIGHT>::getValues(D, z, np, alpha, beta); break;
        case POLYTYPE_GAUSS_LOBATTO:     Polylib::Serial::Derivative<POLYTYPE_GAUSS_LOBATTO>    ::getValues(D, z, np, alpha, beta); break;
        default:
          INTREPID2_TEST_FOR_ABORT(true,
                                   ">>> ERROR (Polylib::Serial::getDerivative): Not supported poly type.");
          break;
        }
      }

      // -----------------------------------------------------------------------
      // Lagrangian Interpolants 
      // -----------------------------------------------------------------------

      template<EPolyType polyType>
      struct LagrangianInterpolant {
        template<typename zViewType>
        KOKKOS_INLINE_FUNCTION
        static typename zViewType::value_type
        getValue(const ordinal_type i, 
                 const typename zViewType::value_type z, 
                 const zViewType zgj,
                 const ordinal_type np,
                 const double alpha,
                 const double beta);
      };

      template<typename zViewType>
      KOKKOS_INLINE_FUNCTION
      static typename zViewType::value_type
      getLagrangianInterpolant(const ordinal_type i, 
                               const typename zViewType::value_type z, 
                               const zViewType zgj,
                               const ordinal_type np,
                               const double alpha,
                               const double beta,
                               const EPolyType poly) {
        typename zViewType::value_type r_val = 0;
        switch (poly) {
        case POLYTYPE_GAUSS:             r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS>            ::getValue(i, z, zgj, np, alpha, beta); break;
        case POLYTYPE_GAUSS_RADAU_LEFT:  r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_RADAU_LEFT> ::getValue(i, z, zgj, np, alpha, beta); break;
        case POLYTYPE_GAUSS_RADAU_RIGHT: r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_RADAU_RIGHT>::getValue(i, z, zgj, np, alpha, beta); break;
        case POLYTYPE_GAUSS_LOBATTO:     r_val = Polylib::Serial::LagrangianInterpolant<POLYTYPE_GAUSS_LOBATTO>    ::getValue(i, z, zgj, np, alpha, beta); break;
        default:
          INTREPID2_TEST_FOR_ABORT(true,
                                   ">>> ERROR (Polylib::Serial::getLagrangianInterpolant): Not supported poly type.");
          break;
        }
        return r_val;
      }

      // -----------------------------------------------------------------------
      // Interpolation operators 
      // -----------------------------------------------------------------------

      template<EPolyType polyType>
      struct InterpolationOperator {
        template<typename imViewType,
                 typename zgrjViewType,
                 typename zmViewType>
        KOKKOS_INLINE_FUNCTION
        static void 
        getMatrix(      imViewType im,
                  const zgrjViewType zgrj,
                  const zmViewType zm,
                  const ordinal_type nz,
                  const ordinal_type mz,
                  const double alpha,
                  const double beta);
      };

      template<typename imViewType,
               typename zgrjViewType,
               typename zmViewType>
      KOKKOS_INLINE_FUNCTION
      static void
      getInterpolationOperator(      imViewType im,
                               const zgrjViewType zgrj,
                               const zmViewType zm,
                               const ordinal_type nz,
                               const ordinal_type mz,
                               const double alpha,
                               const double beta,
                               const EPolyType poly) {
        switch (poly) {
        case POLYTYPE_GAUSS:             Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS>            ::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
        case POLYTYPE_GAUSS_RADAU_LEFT:  Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS_RADAU_LEFT> ::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
        case POLYTYPE_GAUSS_RADAU_RIGHT: Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS_RADAU_RIGHT>::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
        case POLYTYPE_GAUSS_LOBATTO:     Polylib::Serial::InterpolationOperator<POLYTYPE_GAUSS_LOBATTO>    ::getMatrix(im, zgrj, zm, nz, mz, alpha, beta); break;
        default:
          INTREPID2_TEST_FOR_ABORT(true,
                                   ">>> ERROR (Polylib::Serial::getInterpolationOperator): Not supported poly type.");
          break;
        }
      }

      
      // -----------------------------------------------------------------------
      // Polynomial functions 
      // -----------------------------------------------------------------------

      template<typename zViewType,
               typename polyiViewType,
               typename polydViewType>
      KOKKOS_INLINE_FUNCTION
      static void 
      JacobiPolynomial(const ordinal_type np, 
                       const zViewType z, 
                             polyiViewType poly_in, 
                             polydViewType polyd,
                       const ordinal_type n, 
                       const double alpha, 
                       const double beta);


      template<typename zViewType,
               typename polydViewType>
      KOKKOS_INLINE_FUNCTION
      static void 
      JacobiPolynomialDerivative(const ordinal_type np, 
                                 const zViewType z, 
                                       polydViewType polyd, 
                                 const ordinal_type n,
                                 const double alpha, 
                                 const double beta);

      // -----------------------------------------------------------------------
      // Helper functions. 
      // -----------------------------------------------------------------------
      
      template<typename zViewType,
               bool DeflationEnabled = false>
      KOKKOS_INLINE_FUNCTION
      static void   
      JacobiZeros (      zViewType z,
                   const ordinal_type n, 
                   const double alpha, 
                   const double beta);

      template<typename zViewType>
      KOKKOS_INLINE_FUNCTION
      static void
      JacobiZerosPolyDeflation(      zViewType z,
                               const ordinal_type n,
                               const double alpha,
                               const double beta);
      
      template<typename aViewType>
      KOKKOS_INLINE_FUNCTION
      static void
      JacobiZerosTriDiagonal(      aViewType a,
                             const ordinal_type n,
                             const double alpha,
                             const double beta); 

      template<typename aViewType,
               typename bViewType>
      KOKKOS_INLINE_FUNCTION
      static void   
      JacobiZeros(aViewType a,
                  bViewType b,
                  const ordinal_type n, 
                  const double alpha, 
                  const double beta);


      template<typename dViewType,
               typename eViewType>
      KOKKOS_INLINE_FUNCTION      
      static void   
      TriQL(      dViewType d, 
                  eViewType e,
            const ordinal_type n);


      KOKKOS_INLINE_FUNCTION      
      static double
      GammaFunction(const double x);

    };

    // -----------------------------------------------------------------------
  };
  
} // end of Intrepid namespace

  // include templated definitions
#include <Intrepid2_PolylibDef.hpp>

#endif