47#ifndef __INTREPID2_PROJECTIONTOOLS_HPP__
48#define __INTREPID2_PROJECTIONTOOLS_HPP__
50#include "Intrepid2_ConfigDefs.hpp"
54#include "Shards_CellTopology.hpp"
55#include "Shards_BasicTopologies.hpp"
100#include "Teuchos_LAPACK.hpp"
105#ifdef HAVE_INTREPID2_KOKKOSKERNELS
106#include "KokkosBatched_QR_Serial_Internal.hpp"
107#include "KokkosBatched_ApplyQ_Serial_Internal.hpp"
108#include "KokkosBatched_Trsv_Serial_Internal.hpp"
109#include "KokkosBatched_Util.hpp"
114namespace Experimental {
182template<
typename DeviceType>
185 using ExecSpaceType =
typename DeviceType::execution_space;
186 using MemSpaceType =
typename DeviceType::memory_space;
187 using EvalPointsType =
typename ProjectionStruct<DeviceType, double>::EvalPointsType;
206 template<
typename BasisType,
207 typename ortValueType,
class ...ortProperties>
210 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
211 const BasisType* cellBasis,
213 const EvalPointsType evalPointType = EvalPointsType::TARGET
235 template<
typename basisCoeffsValueType,
class ...basisCoeffsProperties,
236 typename funValsValueType,
class ...funValsProperties,
238 typename ortValueType,
class ...ortProperties>
240 getL2BasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
241 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
242 const typename BasisType::ScalarViewType evaluationPoints,
243 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
244 const BasisType* cellBasis,
264 template<
typename BasisType>
267 const BasisType* cellBasis,
269 const EvalPointsType evalPointType = EvalPointsType::TARGET
295 template<
typename basisCoeffsValueType,
class ...basisCoeffsProperties,
296 typename funValsValueType,
class ...funValsProperties,
298 typename ortValueType,
class ...ortProperties>
300 getL2DGBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
301 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
302 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
303 const BasisType* cellBasis,
328 template<
typename basisViewType,
typename targetViewType,
typename BasisType>
331 const targetViewType targetAtTargetEPoints,
332 const BasisType* cellBasis,
355 template<
typename BasisType,
typename OrientationViewType >
358 typename BasisType::ScalarViewType gradEvalPoints,
359 const OrientationViewType cellOrientations,
360 const BasisType* cellBasis,
362 const EvalPointsType evalPointType = EvalPointsType::TARGET
388 template<
class BasisCoeffsViewType,
class TargetValueViewType,
class TargetGradViewType,
389 class BasisType,
class OrientationViewType>
392 const TargetValueViewType targetAtEvalPoints,
393 const TargetGradViewType targetGradAtGradEvalPoints,
394 const typename BasisType::ScalarViewType evaluationPoints,
395 const typename BasisType::ScalarViewType gradEvalPoints,
396 const OrientationViewType cellOrientations,
397 const BasisType* cellBasis,
420 template<
typename BasisType,
421 typename ortValueType,
class ...ortProperties>
424 typename BasisType::ScalarViewType curlEvalPoints,
425 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
426 const BasisType* cellBasis,
428 const EvalPointsType evalPointType = EvalPointsType::TARGET
456 template<
typename basisCoeffsValueType,
class ...basisCoeffsProperties,
457 typename funValsValueType,
class ...funValsProperties,
459 typename ortValueType,
class ...ortProperties>
461 getHCurlBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
462 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
463 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetCurlAtCurlEvalPoints,
464 const typename BasisType::ScalarViewType evaluationPoints,
465 const typename BasisType::ScalarViewType curlEvalPoints,
466 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
467 const BasisType* cellBasis,
490 template<
typename BasisType,
491 typename ortValueType,
class ...ortProperties>
494 typename BasisType::ScalarViewType divEvalPoints,
495 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
496 const BasisType* cellBasis,
498 const EvalPointsType evalPointType = EvalPointsType::TARGET
524 template<
typename basisCoeffsValueType,
class ...basisCoeffsProperties,
525 typename funValsValueType,
class ...funValsProperties,
527 typename ortValueType,
class ...ortProperties>
529 getHDivBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
530 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
531 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetDivAtDivEvalPoints,
532 const typename BasisType::ScalarViewType evaluationPoints,
533 const typename BasisType::ScalarViewType divEvalPoints,
534 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
535 const BasisType* cellBasis,
554 template<
typename BasisType,
555 typename ortValueType,
class ...ortProperties>
558 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
559 const BasisType* cellBasis,
561 const EvalPointsType evalPointType = EvalPointsType::TARGET
582 template<
typename basisCoeffsValueType,
class ...basisCoeffsProperties,
583 typename funValsValueType,
class ...funValsProperties,
585 typename ortValueType,
class ...ortProperties>
587 getHVolBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
588 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
589 const typename BasisType::ScalarViewType evaluationPoints,
590 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
591 const BasisType* cellBasis,
607 std::string systemName_;
608 bool matrixIndependentOfCell_;
617 ElemSystem (std::string systemName,
bool matrixIndependentOfCell) :
618 systemName_(systemName), matrixIndependentOfCell_(matrixIndependentOfCell){};
647 template<
typename ViewType1,
typename ViewType2,
typename ViewType3,
typename ViewType4>
648 void solve(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 tau,
649 ViewType3 w,
const ViewType4 elemDof, ordinal_type n, ordinal_type m=0) {
650#ifdef HAVE_INTREPID2_KOKKOSKERNELS
651 solveDevice(basisCoeffs, elemMat, elemRhs, tau,
654 solveHost(basisCoeffs, elemMat, elemRhs, tau,
662#ifdef HAVE_INTREPID2_KOKKOSKERNELS
663 template<
typename ViewType1,
typename ViewType2,
typename ViewType3,
typename ViewType4>
664 void solveDevice(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 taul,
665 ViewType3 work,
const ViewType4 elemDof, ordinal_type n, ordinal_type m) {
666 using HostDeviceType = Kokkos::Device<Kokkos::DefaultHostExecutionSpace,Kokkos::HostSpace>;
668 ordinal_type numCells = basisCoeffs.extent(0);
670 if(matrixIndependentOfCell_) {
671 auto A0 = Kokkos::subview(elemMat, 0, Kokkos::ALL(), Kokkos::ALL());
672 auto tau0 = Kokkos::subview(taul, 0, Kokkos::ALL());
674 Kokkos::DynRankView<typename ViewType2::value_type, HostDeviceType> A0_host(
"A0_host", A0.extent(0),A0.extent(1));
675 auto A0_device = Kokkos::create_mirror_view(
typename DeviceType::memory_space(), A0_host);
676 Kokkos::deep_copy(A0_device, A0);
677 Kokkos::deep_copy(A0_host, A0_device);
679 for(ordinal_type i=n; i<n+m; ++i)
680 for(ordinal_type j=0; j<n; ++j)
681 A0_host(i,j) = A0_host(j,i);
683 Kokkos::DynRankView<typename ViewType2::value_type, HostDeviceType> tau0_host(
"A0_host", tau0.extent(0));
684 auto tau0_device = Kokkos::create_mirror_view(
typename DeviceType::memory_space(), tau0_host);
685 auto w0_host = Kokkos::create_mirror_view(Kokkos::subview(work, 0, Kokkos::ALL()));
688 KokkosBatched::SerialQR_Internal::invoke(A0_host.extent(0), A0_host.extent(1),
689 A0_host.data(), A0_host.stride_0(), A0_host.stride_1(),
690 tau0_host.data(), tau0_host.stride_0(), w0_host.data());
692 Kokkos::deep_copy(A0_device, A0_host);
693 Kokkos::deep_copy(A0, A0_device);
694 Kokkos::deep_copy(tau0_device, tau0_host);
695 Kokkos::deep_copy(tau0, tau0_device);
697 Kokkos::parallel_for (systemName_,
698 Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
699 KOKKOS_LAMBDA (
const size_t ic) {
700 auto w = Kokkos::subview(work, ic, Kokkos::ALL());
702 auto b = Kokkos::subview(elemRhs, ic, Kokkos::ALL());
705 KokkosBatched::SerialApplyQ_RightForwardInternal::invoke(
706 1, A0.extent(0), A0.extent(1),
707 A0.data(), A0.stride_0(), A0.stride_1(),
708 tau0.data(), tau0.stride_0(),
709 b.data(), 1, b.stride_0(),
713 KokkosBatched::SerialTrsvInternalUpper<KokkosBatched::Algo::Trsv::Unblocked>::invoke(
false,
716 A0.data(), A0.stride_0(), A0.stride_1(),
717 b.data(), b.stride_0());
720 for(ordinal_type i=0; i<n; ++i){
721 basisCoeffs(ic,elemDof(i)) = b(i);
727 Kokkos::parallel_for (systemName_,
728 Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
729 KOKKOS_LAMBDA (
const size_t ic) {
731 auto A = Kokkos::subview(elemMat, ic, Kokkos::ALL(), Kokkos::ALL());
732 auto tau = Kokkos::subview(taul, ic, Kokkos::ALL());
733 auto w = Kokkos::subview(work, ic, Kokkos::ALL());
735 for(ordinal_type i=n; i<n+m; ++i)
736 for(ordinal_type j=0; j<n; ++j)
740 KokkosBatched::SerialQR_Internal::invoke(A.extent(0), A.extent(1),
741 A.data(), A.stride_0(), A.stride_1(), tau.data(), tau.stride_0(), w.data());
743 auto b = Kokkos::subview(elemRhs, ic, Kokkos::ALL());
746 KokkosBatched::SerialApplyQ_RightForwardInternal::invoke(
747 1, A.extent(0), A.extent(1),
748 A.data(), A.stride_0(), A.stride_1(),
749 tau.data(), tau.stride_0(),
750 b.data(), 1, b.stride_0(),
754 KokkosBatched::SerialTrsvInternalUpper<KokkosBatched::Algo::Trsv::Unblocked>::invoke(
false,
757 A.data(), A.stride_0(), A.stride_1(),
758 b.data(), b.stride_0());
761 for(ordinal_type i=0; i<n; ++i){
762 basisCoeffs(ic,elemDof(i)) = b(i);
772 template<
typename ViewType1,
typename ViewType2,
typename ViewType3,
typename ViewType4>
773 void solveHost(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 ,
774 ViewType3,
const ViewType4 elemDof, ordinal_type n, ordinal_type m) {
775 using value_type =
typename ViewType2::value_type;
776 using device_type = DeviceType;
777 using host_exec_space = Kokkos::DefaultHostExecutionSpace;
778 using host_memory_space = Kokkos::HostSpace;
779 using host_device_type = Kokkos::Device<host_exec_space,host_memory_space>;
780 using vector_host_type = Kokkos::View<value_type*,host_device_type>;
781 using scratch_host_type = Kokkos::View<value_type*,host_exec_space::scratch_memory_space>;
782 using matrix_host_type = Kokkos::View<value_type**,Kokkos::LayoutLeft,host_device_type>;
783 using do_not_init_tag = Kokkos::ViewAllocateWithoutInitializing;
784 using host_team_policy_type = Kokkos::TeamPolicy<host_exec_space>;
785 using host_range_policy_type = Kokkos::RangePolicy<host_exec_space>;
791 const ordinal_type numCells = basisCoeffs.extent(0);
792 const ordinal_type numRows = m+n, numCols = n;
795 Teuchos::LAPACK<ordinal_type,value_type> lapack;
798 Kokkos::View<ordinal_type*,host_device_type> elemDof_host(do_not_init_tag(
"elemDof_host"), elemDof.extent(0));
800 auto elemDof_device = Kokkos::create_mirror_view(
typename device_type::memory_space(), elemDof_host);
801 Kokkos::deep_copy(elemDof_device, elemDof); Kokkos::fence();
802 Kokkos::deep_copy(elemDof_host, elemDof_device);
806 auto elemRhs_host = Kokkos::create_mirror_view_and_copy(host_memory_space(), elemRhs);
807 auto elemMat_host = Kokkos::create_mirror_view_and_copy(host_memory_space(), elemMat);
810 auto basisCoeffs_host = Kokkos::create_mirror_view_and_copy(host_memory_space(), basisCoeffs);
812 if (matrixIndependentOfCell_) {
814 matrix_host_type A(do_not_init_tag(
"A"), numRows, numRows);
816 for (ordinal_type j=0;j<numRows;++j)
817 for (ordinal_type i=0;i<numRows;++i)
818 A(i, j) = elemMat_host(0, i, j);
820 for (ordinal_type j=0;j<numCols;++j)
821 for (ordinal_type i=numCols;i<numRows;++i)
825 ordinal_type lwork(-1);
827 ordinal_type info(0);
830 numRows, numRows, numCells,
838 matrix_host_type C(do_not_init_tag(
"C"), numRows, numCells);
840 host_range_policy_type policy(0, numCells);
843 (
"ProjectionTools::solveHost::matrixIndependentOfCell::pack",
844 policy, [=](
const ordinal_type & ic) {
845 for (ordinal_type i=0;i<numRows;++i)
846 C(i,ic) = elemRhs_host(ic, i);
851 vector_host_type work(do_not_init_tag(
"work"), lwork);
852 ordinal_type info(0);
854 numRows, numRows, numCells,
859 INTREPID2_TEST_FOR_EXCEPTION
860 (info != 0, std::runtime_error,
"GELS return non-zero info code");
864 (
"ProjectionTools::solveHost::matrixIndependentOfCell::unpack",
865 policy, [=](
const ordinal_type & ic) {
866 for (ordinal_type i=0;i<numCols;++i)
867 basisCoeffs_host(ic,elemDof_host(i)) = C(i,ic);
871 const ordinal_type level(0);
872 host_team_policy_type policy(numCells, 1, 1);
875 ordinal_type lwork(-1);
877 ordinal_type info(0);
888 const ordinal_type per_team_extent = numRows*numRows + numRows + lwork;
889 const ordinal_type per_team_scratch = scratch_host_type::shmem_size(per_team_extent);
890 policy.set_scratch_size(level, Kokkos::PerTeam(per_team_scratch));
894 (
"ProjectionTools::solveHost::matrixDependentOfCell",
895 policy, [=](
const typename host_team_policy_type::member_type& member) {
896 const ordinal_type ic = member.league_rank();
898 scratch_host_type scratch(member.team_scratch(level), per_team_extent);
899 value_type * sptr = scratch.data();
902 matrix_host_type A(sptr, numRows, numRows); sptr += A.span();
903 for (ordinal_type j=0;j<numRows;++j)
904 for (ordinal_type i=0;i<numRows;++i)
905 A(i, j) = elemMat_host(ic, i, j);
907 for (ordinal_type j=0;j<numCols;++j)
908 for (ordinal_type i=numCols;i<numRows;++i)
911 vector_host_type c(sptr, numRows); sptr += c.span();
912 for (ordinal_type i=0;i<numRows;++i)
913 c(i) = elemRhs_host(ic, i);
915 vector_host_type work(sptr, lwork); sptr += work.span();
916 ordinal_type info(0);
923 INTREPID2_TEST_FOR_EXCEPTION
924 (info != 0, std::runtime_error,
"GELS return non-zero info code");
927 for (ordinal_type i=0;i<numCols;++i)
928 basisCoeffs_host(ic,elemDof_host(i)) = c(i);
931 Kokkos::deep_copy(basisCoeffs, basisCoeffs_host);
Header file for the abstract base class Intrepid2::Basis.
Header file for the Intrepid2::Basis_HCURL_HEX_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_HEX_In_FEM class.
Header file for the Intrepid2::Basis_HCURL_QUAD_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_QUAD_In_FEM class.
Header file for the Intrepid2::Basis_HCURL_TET_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_TET_In_FEM class.
Header file for the Intrepid2::Basis_HCURL_TRI_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_TRI_In_FEM class.
Header file for the Intrepid2::Basis_HCURL_WEDGE_I1_FEM class.
Header file for the Intrepid2::Basis_HDIV_HEX_I1_FEM class.
Header file for the Intrepid2::Basis_HDIV_HEX_In_FEM class.
Header file for the Intrepid2::Basis_HDIV_QUAD_I1_FEM class.
Header file for the Intrepid2::Basis_HDIV_QUAD_In_FEM class.
Header file for the Intrepid2::Basis_HDIV_TET_I1_FEM class.
Header file for the Intrepid2::Basis_HDIV_TET_In_FEM class.
Header file for the Intrepid2::Basis_HDIV_TRI_I1_FEM class.
Header file for the Intrepid2::Basis_HDIV_TRI_In_FEM class.
Header file for the Intrepid2::Basis_HDIV_WEDGE_I1_FEM class.
Header file for the Intrepid2::Basis_HGRAD_HEX_Cn_FEM class.
Header file for the Intrepid2::Basis_HGRAD_LINE_Cn_FEM class.
Header file for the Intrepid2::Basis_HGRAD_QUAD_Cn_FEM class.
Header file for the Intrepid2::Basis_HGRAD_TET_Cn_FEM class.
Header file for the Intrepid2::Basis_HGRAD_TRI_Cn_FEM class.
Header file for the Intrepid2::Basis_HVOL_LINE_Cn_FEM class.
Header file for the Intrepid2::Basis_HVOL_TRI_Cn_FEM class.
Header file for the Intrepid2::Experimental::ProjectionStruct.
Contains definitions of custom data types in Intrepid2.
Header function for Intrepid2::Util class and other utility functions.
An helper class to compute the evaluation points and weights needed for performing projections.