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" 114 namespace Experimental {
182 template<
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 solveParallel(basisCoeffs, elemMat, elemRhs, tau,
662 #ifdef HAVE_INTREPID2_KOKKOSKERNELS 663 template<
typename ViewType1,
typename ViewType2,
typename ViewType3,
typename ViewType4>
664 void solveParallel(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 taul,
665 ViewType3 work,
const ViewType4 elemDof, ordinal_type n, ordinal_type m) {
666 using HostSpaceType =
typename Kokkos::Impl::is_space<DeviceType>::host_mirror_space::execution_space;
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, HostSpaceType> 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, HostSpaceType> 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 solveSerial(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 ,
774 ViewType3,
const ViewType4 elemDof, ordinal_type n, ordinal_type m) {
775 using valueType =
typename ViewType2::value_type;
776 using HostSpaceType =
typename Kokkos::Impl::is_space<DeviceType>::host_mirror_space::execution_space;
777 Kokkos::View<valueType**,Kokkos::LayoutLeft,HostSpaceType>
778 serialElemMat(
"serialElemMat", n+m, n+m);
779 Teuchos::LAPACK<ordinal_type,valueType> lapack_;
780 ordinal_type numCells = basisCoeffs.extent(0);
782 if(matrixIndependentOfCell_) {
783 ViewType2 elemRhsTrans(
"transRhs", elemRhs.extent(1), elemRhs.extent(0));
784 Kokkos::View<valueType**,Kokkos::LayoutLeft,HostSpaceType>
785 pivVec(
"pivVec", m+n + std::max(m+n, numCells), 1);
787 Kokkos::View<valueType**,Kokkos::LayoutLeft,HostSpaceType> serialElemRhs(
"serialElemRhs", n+m, numCells);
789 Kokkos::DynRankView<typename ViewType2::value_type, HostSpaceType> A_host(
"A0_host", elemMat.extent(1),elemMat.extent(2));
790 auto A_device = Kokkos::create_mirror_view(
typename DeviceType::memory_space(), A_host);
791 Kokkos::deep_copy(A_device, Kokkos::subview(elemMat, 0, Kokkos::ALL(), Kokkos::ALL()));
792 Kokkos::deep_copy(A_host, A_device);
794 auto b = Kokkos::create_mirror_view_and_copy(HostSpaceType(), elemRhs);
796 auto serialBasisCoeffs = Kokkos::create_mirror_view_and_copy(
797 HostSpaceType(), basisCoeffs);
799 for(ordinal_type i=0; i<m+n; ++i) {
800 for(ordinal_type ic=0; ic< numCells; ++ic)
801 serialElemRhs(i,ic) = b(ic,i);
802 for(ordinal_type j=0; j<n; ++j)
803 serialElemMat(j,i) = A_host(j,i);
806 for(ordinal_type i=n; i<n+m; ++i)
807 for(ordinal_type j=0; j<n; ++j)
808 serialElemMat(i,j) = serialElemMat(j,i);
810 ordinal_type info = 0;
811 lapack_.GELS(
'N', n+m, n+m, numCells,
812 serialElemMat.data(), serialElemMat.stride_1(),
813 serialElemRhs.data(), serialElemRhs.stride_1(),
814 pivVec.data(), pivVec.extent(0),
817 for(ordinal_type i=0; i<n; ++i) {
818 for (ordinal_type ic = 0; ic < numCells; ic++)
819 serialBasisCoeffs(ic,elemDof(i)) = serialElemRhs(i,ic);
823 Kokkos::View<valueType**,Kokkos::LayoutLeft,HostSpaceType> pivVec(
"pivVec", 2*(m+n), 1);
824 Kokkos::View<valueType**,Kokkos::LayoutLeft,HostSpaceType> serialElemRhs(
"serialElemRhs", n+m, 1 );
825 for (ordinal_type ic = 0; ic < numCells; ic++) {
826 auto A = Kokkos::create_mirror_view_and_copy(HostSpaceType(),
827 Kokkos::subview(elemMat, ic, Kokkos::ALL(), Kokkos::ALL()));
828 auto b = Kokkos::create_mirror_view_and_copy(HostSpaceType(),
829 Kokkos::subview(elemRhs, ic, Kokkos::ALL()));
830 auto basisCoeffs_ = Kokkos::subview(basisCoeffs, ic, Kokkos::ALL());
831 auto serialBasisCoeffs = Kokkos::create_mirror_view_and_copy(HostSpaceType(),
834 Kokkos::deep_copy(serialElemMat,valueType(0));
836 for(ordinal_type i=0; i<m+n; ++i) {
837 serialElemRhs(i,0) = b(i);
838 for(ordinal_type j=0; j<n; ++j)
839 serialElemMat(j,i) = A(j,i);
842 for(ordinal_type i=n; i<n+m; ++i)
843 for(ordinal_type j=0; j<n; ++j)
844 serialElemMat(i,j) = serialElemMat(j,i);
847 ordinal_type info = 0;
848 lapack_.GELS(
'N', n+m, n+m, 1,
849 serialElemMat.data(), serialElemMat.stride_1(),
850 serialElemRhs.data(), serialElemRhs.stride_1(),
851 pivVec.data(), pivVec.extent(0),
855 std::stringstream ss;
856 ss <<
">>> ERROR (Intrepid::ProjectionTools::getBasisCoeffs): " 857 <<
"LAPACK return with error code: " 859 INTREPID2_TEST_FOR_EXCEPTION(
true, std::runtime_error, ss.str().c_str() );
862 for(ordinal_type i=0; i<n; ++i) {
863 serialBasisCoeffs(elemDof(i)) = serialElemRhs(i,0);
865 Kokkos::deep_copy(basisCoeffs_,serialBasisCoeffs);
Header file for the Intrepid2::Basis_HGRAD_LINE_Cn_FEM class.
Header file for the Intrepid2::Basis_HDIV_TET_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_HEX_In_FEM class.
Header file for the Intrepid2::Basis_HDIV_HEX_I1_FEM class.
Header file for the Intrepid2::Basis_HDIV_WEDGE_I1_FEM class.
Header file for the Intrepid2::Basis_HDIV_HEX_In_FEM class.
Header function for Intrepid2::Util class and other utility functions.
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_HDIV_QUAD_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_TET_In_FEM class.
Header file for the Intrepid2::Basis_HGRAD_TET_Cn_FEM class.
Header file for the Intrepid2::Basis_HDIV_TRI_I1_FEM class.
Header file for the Intrepid2::Basis_HDIV_QUAD_In_FEM class.
Header file for the Intrepid2::Basis_HCURL_QUAD_In_FEM class.
Header file for the Intrepid2::Basis_HDIV_TRI_In_FEM class.
Contains definitions of custom data types in Intrepid2.
Header file for the Intrepid2::Basis_HGRAD_QUAD_Cn_FEM class.
Header file for the Intrepid2::Basis_HDIV_TET_In_FEM class.
Header file for the Intrepid2::Experimental::ProjectionStruct.
Header file for the Intrepid2::Basis_HCURL_QUAD_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_WEDGE_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_TET_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_HEX_I1_FEM class.
An helper class to compute the evaluation points and weights needed for performing projections...
Header file for the Intrepid2::Basis_HCURL_TRI_I1_FEM class.
Header file for the Intrepid2::Basis_HVOL_TRI_Cn_FEM class.
Header file for the abstract base class Intrepid2::Basis.
Header file for the Intrepid2::Basis_HCURL_TRI_In_FEM class.
Header file for the Intrepid2::Basis_HGRAD_HEX_Cn_FEM class.