Intrepid2_ProjectionTools.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

#ifndef __INTREPID2_PROJECTIONTOOLS_HPP__
#define __INTREPID2_PROJECTIONTOOLS_HPP__

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

#include "Shards_CellTopology.hpp"
#include "Shards_BasicTopologies.hpp"

#include "Intrepid2_PointTools.hpp"

#include "Intrepid2_Basis.hpp"

// -- HGRAD family
#include "Intrepid2_HGRAD_LINE_Cn_FEM.hpp"
#include "Intrepid2_HGRAD_QUAD_Cn_FEM.hpp"
#include "Intrepid2_HGRAD_HEX_Cn_FEM.hpp"

#include "Intrepid2_HGRAD_TRI_Cn_FEM.hpp"
#include "Intrepid2_HGRAD_TET_Cn_FEM.hpp"

// -- HCURL family
#include "Intrepid2_HCURL_QUAD_In_FEM.hpp"
#include "Intrepid2_HCURL_HEX_In_FEM.hpp"

#include "Intrepid2_HCURL_TRI_In_FEM.hpp"
#include "Intrepid2_HCURL_TET_In_FEM.hpp"
#include "Intrepid2_HVOL_LINE_Cn_FEM.hpp"

// -- HDIV family
#include "Intrepid2_HDIV_QUAD_In_FEM.hpp"
#include "Intrepid2_HDIV_HEX_In_FEM.hpp"

#include "Intrepid2_HDIV_TRI_In_FEM.hpp"
#include "Intrepid2_HDIV_TET_In_FEM.hpp"
#include "Intrepid2_HVOL_TRI_Cn_FEM.hpp"

// -- Lower order family
#include "Intrepid2_HCURL_QUAD_I1_FEM.hpp"
#include "Intrepid2_HCURL_TRI_I1_FEM.hpp"

#include "Intrepid2_HDIV_QUAD_I1_FEM.hpp"
#include "Intrepid2_HDIV_TRI_I1_FEM.hpp"

#include "Intrepid2_HCURL_HEX_I1_FEM.hpp"
#include "Intrepid2_HCURL_TET_I1_FEM.hpp"
#include "Intrepid2_HCURL_WEDGE_I1_FEM.hpp"

#include "Intrepid2_HDIV_HEX_I1_FEM.hpp"
#include "Intrepid2_HDIV_TET_I1_FEM.hpp"
#include "Intrepid2_HDIV_WEDGE_I1_FEM.hpp"

#include "Teuchos_LAPACK.hpp"
#include "Intrepid2_OrientationTools.hpp"

#include "Intrepid2_ProjectionStruct.hpp"

#ifdef HAVE_INTREPID2_KOKKOSKERNELS
#include "KokkosBatched_QR_Serial_Internal.hpp"
#include "KokkosBatched_ApplyQ_Serial_Internal.hpp"
#include "KokkosBatched_Trsv_Serial_Internal.hpp"
#include "KokkosBatched_Util.hpp"
#endif

namespace Intrepid2 {

namespace Experimental {



template<typename ExecSpaceType>
class ProjectionTools {
public:

  using EvalPointsType = typename ProjectionStruct<ExecSpaceType, double>::EvalPointsType;


  template<typename BasisType,
  typename ortValueType,       class ...ortProperties>
  static void
  getL2EvaluationPoints(typename BasisType::ScalarViewType evaluationPoints,
      const Kokkos::DynRankView<ortValueType,   ortProperties...>  cellOrientations,
      const BasisType* cellBasis,
      ProjectionStruct<ExecSpaceType, typename BasisType::scalarType> * projStruct,
      const EvalPointsType evalPointType = EvalPointsType::TARGET
  );

  template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
  typename funValsValueType, class ...funValsProperties,
  typename BasisType,
  typename ortValueType,       class ...ortProperties>
  static void
  getL2BasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
      const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
      const typename BasisType::ScalarViewType evaluationPoints,
      const Kokkos::DynRankView<ortValueType,   ortProperties...>  cellOrientations,
      const BasisType* cellBasis,
      ProjectionStruct<ExecSpaceType, typename BasisType::scalarType> * projStruct);


  template<typename BasisType,
  typename ortValueType,       class ...ortProperties>
  static void
  getHGradEvaluationPoints(typename BasisType::ScalarViewType evaluationPoints,
      typename BasisType::ScalarViewType gradEvalPoints,
      const Kokkos::DynRankView<ortValueType,   ortProperties...>  cellOrientations,
      const BasisType* cellBasis,
      ProjectionStruct<ExecSpaceType, typename BasisType::scalarType> * projStruct,
      const EvalPointsType evalPointType = EvalPointsType::TARGET
  );

  template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
  typename funValsValueType, class ...funValsProperties,
  typename BasisType,
  typename ortValueType,       class ...ortProperties>
  static void
  getHGradBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
      const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
      const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetGradAtGradEvalPoints,
      const typename BasisType::ScalarViewType evaluationPoints,
      const typename BasisType::ScalarViewType gradEvalPoints,
      const Kokkos::DynRankView<ortValueType,   ortProperties...>  cellOrientations,
      const BasisType* cellBasis,
      ProjectionStruct<ExecSpaceType, typename BasisType::scalarType> * projStruct);


  template<typename BasisType,
  typename ortValueType,       class ...ortProperties>
  static void
  getHCurlEvaluationPoints(typename BasisType::ScalarViewType evaluationPoints,
      typename BasisType::ScalarViewType curlEvalPoints,
      const Kokkos::DynRankView<ortValueType,   ortProperties...>  cellOrientations,
      const BasisType* cellBasis,
      ProjectionStruct<ExecSpaceType, typename BasisType::scalarType> * projStruct,
      const EvalPointsType evalPointType = EvalPointsType::TARGET
  );

  template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
  typename funValsValueType, class ...funValsProperties,
  typename BasisType,
  typename ortValueType,       class ...ortProperties>
  static void
  getHCurlBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
      const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
      const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetCurlAtCurlEvalPoints,
      const typename BasisType::ScalarViewType evaluationPoints,
      const typename BasisType::ScalarViewType curlEvalPoints,
      const Kokkos::DynRankView<ortValueType,   ortProperties...>  cellOrientations,
      const BasisType* cellBasis,
      ProjectionStruct<ExecSpaceType, typename BasisType::scalarType> * projStruct);


  template<typename BasisType,
  typename ortValueType,       class ...ortProperties>
  static void
  getHDivEvaluationPoints(typename BasisType::ScalarViewType evaluationPoints,
      typename BasisType::ScalarViewType divEvalPoints,
      const Kokkos::DynRankView<ortValueType,   ortProperties...>  cellOrientations,
      const BasisType* cellBasis,
      ProjectionStruct<ExecSpaceType, typename BasisType::scalarType> * projStruct,
      const EvalPointsType evalPointType = EvalPointsType::TARGET
  );

  template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
  typename funValsValueType, class ...funValsProperties,
  typename BasisType,
  typename ortValueType,       class ...ortProperties>
  static void
  getHDivBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
      const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
      const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetDivAtDivEvalPoints,
      const typename BasisType::ScalarViewType evaluationPoints,
      const typename BasisType::ScalarViewType divEvalPoints,
      const Kokkos::DynRankView<ortValueType,   ortProperties...>  cellOrientations,
      const BasisType* cellBasis,
      ProjectionStruct<ExecSpaceType, typename BasisType::scalarType> * projStruct);

  template<typename BasisType,
  typename ortValueType,       class ...ortProperties>
  static void
  getHVolEvaluationPoints(typename BasisType::ScalarViewType evaluationPoints,
      const Kokkos::DynRankView<ortValueType, ortProperties...>  cellOrientations,
      const BasisType* cellBasis,
      ProjectionStruct<ExecSpaceType, typename BasisType::scalarType> * projStruct,
      const EvalPointsType evalPointType = EvalPointsType::TARGET
  );

  template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
  typename funValsValueType, class ...funValsProperties,
  typename BasisType,
  typename ortValueType,       class ...ortProperties>
  static void
  getHVolBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
      const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
      const typename BasisType::ScalarViewType evaluationPoints,
      const Kokkos::DynRankView<ortValueType,   ortProperties...>  cellOrientations,
      const BasisType* cellBasis,
      ProjectionStruct<ExecSpaceType, typename BasisType::scalarType> * projStruct);


  struct ElemSystem {


    std::string systemName_;
    bool matrixIndependentOfCell_;

    ElemSystem (std::string systemName, bool matrixIndependentOfCell) :
      systemName_(systemName), matrixIndependentOfCell_(matrixIndependentOfCell){};



    template<typename ViewType1, typename ViewType2, typename ViewType3, typename ViewType4>
    void solve(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 tau,
        ViewType3 w,const  ViewType4 elemDof, ordinal_type n, ordinal_type m=0) {
#ifdef HAVE_INTREPID2_KOKKOSKERNELS
      solveParallel(basisCoeffs, elemMat, elemRhs, tau,
          w, elemDof, n, m);
#else
      solveSerial(basisCoeffs, elemMat, elemRhs, tau,
          w, elemDof, n, m);
#endif

    }

#ifdef HAVE_INTREPID2_KOKKOSKERNELS
    template<typename ViewType1, typename ViewType2, typename ViewType3, typename ViewType4>
    void solveParallel(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 taul,
        ViewType3 work,const  ViewType4 elemDof, ordinal_type n, ordinal_type m) {

      ordinal_type numCells = basisCoeffs.extent(0);

      if(matrixIndependentOfCell_) {
        auto A0 = Kokkos::subview(elemMat, 0, Kokkos::ALL(), Kokkos::ALL());
        auto tau0 = Kokkos::subview(taul, 0, Kokkos::ALL());

        auto A0_host = Kokkos::create_mirror_view_and_copy(typename ExecSpaceType::memory_space(), A0);
        auto tau0_host = Kokkos::create_mirror_view(typename ExecSpaceType::memory_space(), tau0);


        for(ordinal_type i=n; i<n+m; ++i)
          for(ordinal_type j=0; j<n; ++j)
            A0_host(i,j) = A0_host(j,i);

        auto w0_host = Kokkos::create_mirror_view(Kokkos::subview(work, 0, Kokkos::ALL()));

        //computing QR of A0. QR is saved in A0 and tau0
        KokkosBatched::SerialQR_Internal::invoke(A0_host.extent(0), A0_host.extent(1),
            A0_host.data(), A0_host.stride_0(), A0_host.stride_1(),
            tau0_host.data(), tau0_host.stride_0(), w0_host.data());

        Kokkos::deep_copy(A0, A0_host);
        Kokkos::deep_copy(tau0, tau0_host);

        Kokkos::parallel_for (systemName_,
            Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
            KOKKOS_LAMBDA (const size_t ic) {
          auto w = Kokkos::subview(work, ic, Kokkos::ALL());

          auto b = Kokkos::subview(elemRhs, ic, Kokkos::ALL());

          //b'*Q0 -> b
          KokkosBatched::SerialApplyQ_RightNoTransForwardInternal::invoke(
              1, A0.extent(0), A0.extent(1),
              A0.data(),  A0.stride_0(), A0.stride_1(),
              tau0.data(), tau0.stride_0(),
              b.data(),  1, b.stride_0(),
              w.data());

          // R0^{-1} b -> b
          KokkosBatched::SerialTrsvInternalUpper<KokkosBatched::Algo::Trsv::Unblocked>::invoke(false,
              A0.extent(0),
              1.0,
              A0.data(), A0.stride_0(), A0.stride_1(),
              b.data(),  b.stride_0());

          //scattering b into the basis coefficients
          for(ordinal_type i=0; i<n; ++i){
            basisCoeffs(ic,elemDof(i)) = b(i);
          }
        });

      } else {

        Kokkos::parallel_for (systemName_,
            Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
            KOKKOS_LAMBDA (const size_t ic) {

          auto A = Kokkos::subview(elemMat, ic, Kokkos::ALL(), Kokkos::ALL());
          auto tau = Kokkos::subview(taul, ic, Kokkos::ALL());
          auto w = Kokkos::subview(work, ic, Kokkos::ALL());

          for(ordinal_type i=n; i<n+m; ++i)
            for(ordinal_type j=0; j<n; ++j)
              A(i,j) = A(j,i);

          //computing QR of A. QR is saved in A and tau
          KokkosBatched::SerialQR_Internal::invoke(A.extent(0), A.extent(1),
              A.data(), A.stride_0(), A.stride_1(), tau.data(), tau.stride_0(), w.data());

          auto b = Kokkos::subview(elemRhs, ic, Kokkos::ALL());

          //b'*Q -> b
          KokkosBatched::SerialApplyQ_RightNoTransForwardInternal::invoke(
              1, A.extent(0), A.extent(1),
              A.data(),  A.stride_0(), A.stride_1(),
              tau.data(), tau.stride_0(),
              b.data(),  1, b.stride_0(),
              w.data());

          // R^{-1} b -> b
          KokkosBatched::SerialTrsvInternalUpper<KokkosBatched::Algo::Trsv::Unblocked>::invoke(false,
              A.extent(0),
              1.0,
              A.data(), A.stride_0(), A.stride_1(),
              b.data(),  b.stride_0());

          //scattering b into the basis coefficients
          for(ordinal_type i=0; i<n; ++i){
            basisCoeffs(ic,elemDof(i)) = b(i);
          }
        });
      }
    }
#endif

    template<typename ViewType1, typename ViewType2, typename ViewType3, typename ViewType4>
    void solveSerial(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 ,
        ViewType3, const  ViewType4 elemDof, ordinal_type n, ordinal_type m) {
      using valueType = typename ViewType2::value_type;
      using host_space_type = typename Kokkos::Impl::is_space<ExecSpaceType>::host_mirror_space::execution_space;
      Kokkos::View<valueType**,Kokkos::LayoutLeft,host_space_type>
      serialElemMat("serialElemMat", n+m, n+m);
      Teuchos::LAPACK<ordinal_type,valueType> lapack_;
      ordinal_type numCells = basisCoeffs.extent(0);

      if(matrixIndependentOfCell_) {
        ViewType2 elemRhsTrans("transRhs", elemRhs.extent(1), elemRhs.extent(0));
        Kokkos::View<valueType**,Kokkos::LayoutLeft,host_space_type>
        pivVec("pivVec", m+n + std::max(m+n, numCells), 1);

        Kokkos::View<valueType**,Kokkos::LayoutLeft,host_space_type> serialElemRhs("serialElemRhs", n+m, numCells);

        auto A = Kokkos::create_mirror_view_and_copy(typename ExecSpaceType::memory_space(),
            Kokkos::subview(elemMat, 0, Kokkos::ALL(), Kokkos::ALL()));
        auto b = Kokkos::create_mirror_view_and_copy(typename ExecSpaceType::memory_space(), elemRhs);

        auto serialBasisCoeffs = Kokkos::create_mirror_view_and_copy(
            typename ExecSpaceType::memory_space(), basisCoeffs);

        for(ordinal_type i=0; i<m+n; ++i) {
          for(ordinal_type ic=0; ic< numCells; ++ic)
            serialElemRhs(i,ic) = b(ic,i);
          for(ordinal_type j=0; j<n; ++j)
            serialElemMat(j,i) = A(j,i);
        }

        for(ordinal_type i=n; i<n+m; ++i)
          for(ordinal_type j=0; j<n; ++j)
            serialElemMat(i,j) = serialElemMat(j,i);

        ordinal_type info = 0;
        lapack_.GELS('N', n+m, n+m, numCells,
            serialElemMat.data(), serialElemMat.stride_1(),
            serialElemRhs.data(), serialElemRhs.stride_1(),
            pivVec.data(), pivVec.extent(0),
            &info);

        for(ordinal_type i=0; i<n; ++i) {
          for (ordinal_type ic = 0; ic < numCells; ic++)
            serialBasisCoeffs(ic,elemDof(i)) = serialElemRhs(i,ic);
        }
      }
      else {
        Kokkos::View<valueType**,Kokkos::LayoutLeft,host_space_type> pivVec("pivVec", 2*(m+n), 1);
        Kokkos::View<valueType**,Kokkos::LayoutLeft,host_space_type> serialElemRhs("serialElemRhs", n+m, 1 );
        for (ordinal_type ic = 0; ic < numCells; ic++) {
          auto A = Kokkos::create_mirror_view_and_copy(typename ExecSpaceType::memory_space(),
              Kokkos::subview(elemMat, ic, Kokkos::ALL(), Kokkos::ALL()));
          auto b = Kokkos::create_mirror_view_and_copy(typename ExecSpaceType::memory_space(),
              Kokkos::subview(elemRhs, ic, Kokkos::ALL()));
          auto basisCoeffs_ = Kokkos::subview(basisCoeffs, ic, Kokkos::ALL());
          auto serialBasisCoeffs = Kokkos::create_mirror_view_and_copy(typename ExecSpaceType::memory_space(),
              basisCoeffs_);

          Kokkos::deep_copy(serialElemMat,valueType(0));  //LAPACK might overwrite the matrix

          for(ordinal_type i=0; i<m+n; ++i) {
            serialElemRhs(i,0) = b(i);
            for(ordinal_type j=0; j<n; ++j)
              serialElemMat(j,i) = A(j,i);
          }

          for(ordinal_type i=n; i<n+m; ++i)
            for(ordinal_type j=0; j<n; ++j)
              serialElemMat(i,j) = serialElemMat(j,i);

          // Using GELS because the matrix can be close to singular.
          ordinal_type info = 0;
          lapack_.GELS('N', n+m, n+m, 1,
              serialElemMat.data(), serialElemMat.stride_1(),
              serialElemRhs.data(), serialElemRhs.stride_1(),
              pivVec.data(), pivVec.extent(0),
              &info);

          if (info) {
            std::stringstream ss;
            ss << ">>> ERROR (Intrepid::ProjectionTools::getBasisCoeffs): "
                << "LAPACK return with error code: "
                << info;
            INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
          }

          for(ordinal_type i=0; i<n; ++i) {
            serialBasisCoeffs(elemDof(i)) = serialElemRhs(i,0);
          }
          Kokkos::deep_copy(basisCoeffs_,serialBasisCoeffs);
        }
      }
    }
  };

};

} //Experimental
} //Intrepid2


// include templated function definitions
#include "Intrepid2_ProjectionToolsDefL2.hpp"
#include "Intrepid2_ProjectionToolsDefHGRAD.hpp"
#include "Intrepid2_ProjectionToolsDefHCURL.hpp"
#include "Intrepid2_ProjectionToolsDefHDIV.hpp"
#include "Intrepid2_ProjectionToolsDefHVOL.hpp"

#endif