diff --git a/cmake/std/atdm/waterman/tweaks/CUDA_COMMON_TWEAKS.cmake b/cmake/std/atdm/waterman/tweaks/CUDA_COMMON_TWEAKS.cmake index 81924d8ba35b..e69de29bb2d1 100644 --- a/cmake/std/atdm/waterman/tweaks/CUDA_COMMON_TWEAKS.cmake +++ b/cmake/std/atdm/waterman/tweaks/CUDA_COMMON_TWEAKS.cmake @@ -1,2 +0,0 @@ -# This test runs out of CUDA memory all on its own (#3340) -ATDM_SET_ENABLE(PanzerAdaptersSTK_CurlLaplacianExample-ConvTest-Quad-Order-4_DISABLE ON) diff --git a/packages/panzer/adapters-stk/example/CurlLaplacianExample/CMakeLists.txt b/packages/panzer/adapters-stk/example/CurlLaplacianExample/CMakeLists.txt index dec956d43b92..e8c5d21200d3 100644 --- a/packages/panzer/adapters-stk/example/CurlLaplacianExample/CMakeLists.txt +++ b/packages/panzer/adapters-stk/example/CurlLaplacianExample/CMakeLists.txt @@ -27,23 +27,23 @@ FOREACH( ORDER 1 2 3 4 ) TRIBITS_ADD_ADVANCED_TEST( CurlLaplacianExample-ConvTest-Quad-Order-${ORDER} TEST_0 EXEC CurlLaplacianExample - ARGS --use-epetra --use-twod --cell="Quad" --x-elements=4 --y-elements=4 --z-elements=4 --basis-order=${ORDER} - PASS_REGULAR_EXPRESSION "ALL PASSED: Epetra" + ARGS --use-tpetra --use-twod --cell="Quad" --x-elements=4 --y-elements=4 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" NUM_MPI_PROCS 4 OUTPUT_FILE MPE-ConvTest-Quad-${ORDER}-04 TEST_1 EXEC CurlLaplacianExample - ARGS --use-epetra --use-twod --cell="Quad" --x-elements=8 --y-elements=8 --z-elements=4 --basis-order=${ORDER} - PASS_REGULAR_EXPRESSION "ALL PASSED: Epetra" + ARGS --use-tpetra --use-twod --cell="Quad" --x-elements=8 --y-elements=8 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" NUM_MPI_PROCS 4 OUTPUT_FILE MPE-ConvTest-Quad-${ORDER}-08 TEST_2 EXEC CurlLaplacianExample - ARGS --use-epetra --use-twod --cell="Quad" --x-elements=16 --y-elements=16 --z-elements=4 --basis-order=${ORDER} - PASS_REGULAR_EXPRESSION "ALL PASSED: Epetra" + ARGS --use-tpetra --use-twod --cell="Quad" --x-elements=16 --y-elements=16 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" NUM_MPI_PROCS 4 OUTPUT_FILE MPE-ConvTest-Quad-${ORDER}-16 TEST_3 EXEC CurlLaplacianExample - ARGS --use-epetra --use-twod --cell="Quad" --x-elements=32 --y-elements=32 --z-elements=4 --basis-order=${ORDER} - PASS_REGULAR_EXPRESSION "ALL PASSED: Epetra" + ARGS --use-tpetra --use-twod --cell="Quad" --x-elements=32 --y-elements=32 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" NUM_MPI_PROCS 4 OUTPUT_FILE MPE-ConvTest-Quad-${ORDER}-32 TEST_4 CMND python diff --git a/packages/panzer/adapters-stk/example/CurlLaplacianExample/main.cpp b/packages/panzer/adapters-stk/example/CurlLaplacianExample/main.cpp index 5475685d1deb..b2c466981b5a 100644 --- a/packages/panzer/adapters-stk/example/CurlLaplacianExample/main.cpp +++ b/packages/panzer/adapters-stk/example/CurlLaplacianExample/main.cpp @@ -86,7 +86,7 @@ #include "EpetraExt_RowMatrixOut.h" #include "EpetraExt_VectorOut.h" -#include "BelosBlockGmresSolMgr.hpp" +#include "BelosPseudoBlockGmresSolMgr.hpp" #include "BelosTpetraAdapter.hpp" #include "Example_BCStrategy_Factory.hpp" @@ -169,7 +169,7 @@ int main(int argc,char * argv[]) // Build command line processor //////////////////////////////////////////////////// - bool useTpetra = false; + bool useTpetra = true; bool threeD = false; int x_elements=20,y_elements=20,z_elements=20; std::string celltype = "Quad"; // or "Tri" (2d), Hex or Tet (3d) @@ -638,11 +638,10 @@ void solveTpetraSystem(panzer::LinearObjContainer & container) Teuchos::RCP problem(new ProblemType(tp_container.get_A(), tp_container.get_x(), tp_container.get_f())); TEUCHOS_ASSERT(problem->setProblem()); - typedef Belos::BlockGmresSolMgr SolverType; + typedef Belos::PseudoBlockGmresSolMgr SolverType; Teuchos::ParameterList belosList; - belosList.set( "Flexible Gmres", false ); // Flexible Gmres will be used to solve this problem - belosList.set( "Num Blocks", 1000 ); // Maximum number of blocks in Krylov factorization + belosList.set( "Num Blocks", 3000 ); // Maximum number of blocks in Krylov factorization belosList.set( "Block Size", 1 ); // Blocksize to be used by iterative solver belosList.set( "Maximum Iterations", 5000 ); // Maximum number of iterations allowed belosList.set( "Maximum Restarts", 1 ); // Maximum number of restarts allowed diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/point_basis_values_evaluator.cpp b/packages/panzer/adapters-stk/test/evaluator_tests/point_basis_values_evaluator.cpp index 501e8ca656b7..ad02a85de429 100644 --- a/packages/panzer/adapters-stk/test/evaluator_tests/point_basis_values_evaluator.cpp +++ b/packages/panzer/adapters-stk/test/evaluator_tests/point_basis_values_evaluator.cpp @@ -328,7 +328,7 @@ namespace panzer { fm.evaluateFields(workset); - PHX::MDField + PHX::MDField basis(basis_q1->name()+"_"+point_rule->getName()+"_"+"basis",layout->basis); fm.getFieldData(basis); out << basis << std::endl; @@ -341,7 +341,7 @@ namespace panzer { for(int i=0;i<4;i++) { for(int j=0;j<4;j++) { for(unsigned int k=0;kbasis_scalar.extent(2);k++) { - TEST_FLOATING_EQUALITY(bases->basis_scalar(i,j,k),basis(i,j,k).val(),1e-10); + TEST_FLOATING_EQUALITY(bases->basis_scalar(i,j,k),basis(i,j,k),1e-10); } } } @@ -472,16 +472,17 @@ namespace panzer { fm.evaluateFields(workset); - PHX::MDField + // NOTE: basis values are always double, not Fad. + PHX::MDField basis(basis_edge->name()+"_"+point_rule->getName()+"_basis",layout->basis_grad); - PHX::MDField + PHX::MDField curl_basis(basis_edge->name()+"_"+point_rule->getName()+"_curl_basis",layout->basis); - PHX::MDField + PHX::MDField jac_inv("CubaturePoints (Degree=4,volume)_jac_inv",point_rule->dl_tensor); - fm.getFieldData(basis); - fm.getFieldData(curl_basis); - fm.getFieldData(jac_inv); + fm.getFieldData(basis); + fm.getFieldData(curl_basis); + fm.getFieldData(jac_inv); WorksetDetailsAccessor wda; std::size_t basisIndex = panzer::getBasisIndex(layout->name(), workset, wda); @@ -494,7 +495,7 @@ namespace panzer { for(int j=0;j<4;j++) { for(unsigned int k=0;kcurl_basis_scalar.extent(2);k++) { for(unsigned int d=0;dbasis_vector.extent(3);d++) { - TEST_FLOATING_EQUALITY(bases->basis_vector(i,j,k,d),basis(i,j,k,d).val(),1e-10); + TEST_FLOATING_EQUALITY(bases->basis_vector(i,j,k,d),basis(i,j,k,d),1e-10); } } } @@ -503,7 +504,7 @@ namespace panzer { for(int i=0;i<4;i++) { for(int j=0;j<4;j++) { for(unsigned int k=0;kcurl_basis_scalar.extent(2);k++) { - TEST_FLOATING_EQUALITY(bases->curl_basis_scalar(i,j,k),curl_basis(i,j,k).val(),1e-10); + TEST_FLOATING_EQUALITY(bases->curl_basis_scalar(i,j,k),curl_basis(i,j,k),1e-10); } } } diff --git a/packages/panzer/disc-fe/src/Panzer_BasisValues2.cpp b/packages/panzer/disc-fe/src/Panzer_BasisValues2.cpp index 4830c7a15c83..a18be416bee6 100644 --- a/packages/panzer/disc-fe/src/Panzer_BasisValues2.cpp +++ b/packages/panzer/disc-fe/src/Panzer_BasisValues2.cpp @@ -43,1512 +43,22 @@ #include "PanzerDiscFE_config.hpp" #include "Panzer_Traits.hpp" #include "Panzer_BasisValues2.hpp" - -#include "Panzer_CommonArrayFactories.hpp" -#include "Kokkos_ViewFactory.hpp" - -#include "Intrepid2_Utils.hpp" -#include "Intrepid2_FunctionSpaceTools.hpp" -#include "Intrepid2_Orientation.hpp" -#include "Intrepid2_OrientationTools.hpp" - +#include "Panzer_BasisValues2_impl.hpp" namespace panzer { -template -void panzer::BasisValues2:: -evaluateValues(const PHX::MDField & cub_points, - const PHX::MDField & jac, - const PHX::MDField & jac_det, - const PHX::MDField & jac_inv) -{ - PHX::MDField weighted_measure; - PHX::MDField vertex_coordinates; - build_weighted = false; - evaluateValues(cub_points,jac,jac_det,jac_inv,weighted_measure,vertex_coordinates,false); -} - -template -void panzer::BasisValues2:: -evaluateValues(const PHX::MDField & cub_points, - const PHX::MDField & jac, - const PHX::MDField & jac_det, - const PHX::MDField & jac_inv, - const PHX::MDField & weighted_measure, - const PHX::MDField & vertex_coordinates, - bool use_vertex_coordinates) -{ - MDFieldArrayFactory af("",ddims_,true); - - int num_dim = basis_layout->dimension(); - - // currently this just copies on the basis objects are converted - // in intrepid - evaluateReferenceValues(cub_points,compute_derivatives,use_vertex_coordinates); - - PureBasis::EElementSpace elmtspace = getElementSpace(); - if(elmtspace==PureBasis::CONST || - elmtspace==PureBasis::HGRAD) { - Intrepid2::FunctionSpaceTools:: - HGRADtransformVALUE(basis_scalar.get_view(), - basis_ref_scalar.get_view()); - - if(build_weighted) { - Intrepid2::FunctionSpaceTools:: - multiplyMeasure(weighted_basis_scalar.get_view(), - weighted_measure.get_view(), - basis_scalar.get_view()); - } - } - else if(elmtspace==PureBasis::HCURL) { - Intrepid2::FunctionSpaceTools:: - HCURLtransformVALUE(basis_vector.get_view(), - jac_inv.get_view(), - basis_ref_vector.get_view()); - - if(build_weighted) { - Intrepid2::FunctionSpaceTools:: - multiplyMeasure(weighted_basis_vector.get_view(), - weighted_measure.get_view(), - basis_vector.get_view()); - } - } - else if(elmtspace==PureBasis::HDIV) - { - Intrepid2::FunctionSpaceTools:: - HDIVtransformVALUE(basis_vector.get_view(), - jac.get_view(), - jac_det.get_view(), - basis_ref_vector.get_view()); - - if(build_weighted) { - Intrepid2::FunctionSpaceTools:: - multiplyMeasure(weighted_basis_vector.get_view(), - weighted_measure.get_view(), - basis_vector.get_view()); - } - } - else { TEUCHOS_ASSERT(false); } - - if(elmtspace==PureBasis::HGRAD && compute_derivatives) { - Intrepid2::FunctionSpaceTools:: - HGRADtransformGRAD(grad_basis.get_view(), - jac_inv.get_view(), - grad_basis_ref.get_view()); - - if(build_weighted) { - Intrepid2::FunctionSpaceTools:: - multiplyMeasure(weighted_grad_basis.get_view(), - weighted_measure.get_view(), - grad_basis.get_view()); - } - } - else if(elmtspace==PureBasis::HCURL && num_dim==2 && compute_derivatives) { - Intrepid2::FunctionSpaceTools:: - HDIVtransformDIV(curl_basis_scalar.get_view(), - jac_det.get_view(), // note only volume deformation is needed! - // this relates directly to this being in - // the divergence space in 2D! - curl_basis_ref_scalar.get_view()); - - if(build_weighted) { - Intrepid2::FunctionSpaceTools:: - multiplyMeasure(weighted_curl_basis_scalar.get_view(), - weighted_measure.get_view(), - curl_basis_scalar.get_view()); - } - } - else if(elmtspace==PureBasis::HCURL && num_dim==3 && compute_derivatives) { - Intrepid2::FunctionSpaceTools:: - HCURLtransformCURL(curl_basis_vector.get_view(), - jac.get_view(), - jac_det.get_view(), - curl_basis_ref_vector.get_view()); - - if(build_weighted) { - Intrepid2::FunctionSpaceTools:: - multiplyMeasure(weighted_curl_basis_vector.get_view(), - weighted_measure.get_view(), - curl_basis_vector.get_view()); - } - } - else if(elmtspace==PureBasis::HDIV && compute_derivatives) { - Intrepid2::FunctionSpaceTools:: - HDIVtransformDIV(div_basis.get_view(), - jac_det.get_view(), - div_basis_ref.get_view()); - - if(build_weighted) { - Intrepid2::FunctionSpaceTools:: - multiplyMeasure(weighted_div_basis.get_view(), - weighted_measure.get_view(), - div_basis.get_view()); - } - } - - // If basis supports coordinate values at basis points, then - // compute these values - if(use_vertex_coordinates) { - // Teuchos::RCP > coords - // = Teuchos::rcp_dynamic_cast >(intrepid_basis); - // if (!Teuchos::is_null(coords)) { -/* - ArrayDynamic dyn_basis_coordinates_ref = af.buildArray("basis_coordinates_ref",basis_coordinates_ref.extent(0),basis_coordinates_ref.extent(1)); - coords->getDofCoords(dyn_basis_coordinates_ref); - - // fill in basis coordinates - for (size_type i = 0; i < basis_coordinates_ref.extent(0); ++i) - for (size_type j = 0; j < basis_coordinates_ref.extent(1); ++j) - basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j); -*/ - - Intrepid2::CellTools cell_tools; - cell_tools.mapToPhysicalFrame(basis_coordinates.get_view(), - basis_coordinates_ref.get_view(), - vertex_coordinates.get_view(), - intrepid_basis->getBaseCellTopology()); - } -} - - - - - - - -template -void panzer::BasisValues2:: -evaluateBasisCoordinates(const PHX::MDField & vertex_coordinates) -{ - MDFieldArrayFactory af("",ddims_,true); - - // intrepid_basis->getDofCoords requires DynRankView, but basis_coordinates_ref is more of a static View - // We use an auxiliary 'dyn' array to get around this - using coordsScalarType = typename Intrepid2::Basis::scalarType; - auto dyn_basis_coordinates_ref = af.buildArray("basis_coordinates_ref", - basis_coordinates_ref.extent(0), - basis_coordinates_ref.extent(1)); - intrepid_basis->getDofCoords(dyn_basis_coordinates_ref.get_view()); - - // fill in basis coordinates - for (int i = 0; i < basis_coordinates_ref.extent_int(0); ++i) - for (int j = 0; j < basis_coordinates_ref.extent_int(1); ++j) - basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j); - - Intrepid2::CellTools cell_tools; - cell_tools.mapToPhysicalFrame(basis_coordinates.get_view(), - basis_coordinates_ref.get_view(), - vertex_coordinates.get_view(), - intrepid_basis->getBaseCellTopology()); -} - - - - - -template -void panzer::BasisValues2:: -evaluateValues(const PHX::MDField & cub_points, - const PHX::MDField & jac, - const PHX::MDField & jac_det, - const PHX::MDField & jac_inv, - const PHX::MDField & weighted_measure, - const PHX::MDField & vertex_coordinates, - bool use_vertex_coordinates) -{ - - PureBasis::EElementSpace elmtspace = getElementSpace(); - - if(elmtspace == PureBasis::CONST){ - evaluateValues_Const(cub_points,jac_inv,weighted_measure); - } else if(elmtspace == PureBasis::HGRAD){ - evaluateValues_HGrad(cub_points,jac_inv,weighted_measure); - } else if(elmtspace == PureBasis::HCURL){ - evaluateValues_HCurl(cub_points,jac,jac_det,jac_inv,weighted_measure); - } else if(elmtspace == PureBasis::HDIV){ - evaluateValues_HDiv(cub_points,jac,jac_det,weighted_measure); - } else { - TEUCHOS_TEST_FOR_EXCEPT_MSG(true,"panzer::BasisValues2::evaluateValues : Element space not recognized."); - } - - if(use_vertex_coordinates) { - TEUCHOS_TEST_FOR_EXCEPT_MSG(elmtspace == PureBasis::CONST,"panzer::BasisValues2::evaluateValues : Const basis cannot have basis coordinates."); - evaluateBasisCoordinates(vertex_coordinates); - } - -} - -template -void panzer::BasisValues2:: -evaluateValues_Const(const PHX::MDField & cub_points, - const PHX::MDField & jac_inv, - const PHX::MDField & weighted_measure) -{ - - TEUCHOS_ASSERT(getElementSpace() == PureBasis::CONST); - - typedef Intrepid2::FunctionSpaceTools fst; - MDFieldArrayFactory af("",ddims_,true); - - const panzer::PureBasis & basis = *(basis_layout->getBasis()); - - const int num_points = basis_layout->numPoints(); - const int num_basis = basis.cardinality(); - const int num_dim = basis_layout->dimension(); - const int num_cells = basis_layout->numCells(); - - auto cell_basis_scalar = af.buildStaticArray("cell_basis_scalar",1,num_basis,num_points); - auto cell_cub_points = af.buildStaticArray("cell_cub_points",num_points,num_dim); - auto cell_grad_basis = af.buildStaticArray("cell_grad_basis",1,num_basis,num_points,num_dim); - auto cell_jac_inv = af.buildStaticArray("cell_jac_inv",1,num_points,num_dim,num_dim); - - auto cell_basis_ref_scalar = af.buildStaticArray("cell_basis_ref_scalar",num_basis,num_points); - auto cell_grad_basis_ref = af.buildStaticArray("cell_grad_basis_ref",num_basis,num_points,num_dim); - - for(int cell=0;cellgetValues(cell_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE); - - if(compute_derivatives){ - Kokkos::deep_copy(cell_grad_basis_ref.get_view(),0.0); - } - - // ============================================= - // Transform reference values to physical values - - fst::HGRADtransformVALUE(cell_basis_scalar.get_view(),cell_basis_ref_scalar.get_view()); - for(int b=0;b -void panzer::BasisValues2:: -evaluateValues_HGrad(const PHX::MDField & cub_points, - const PHX::MDField & jac_inv, - const PHX::MDField & weighted_measure) -{ - - TEUCHOS_ASSERT(getElementSpace() == PureBasis::HGRAD); - - typedef Intrepid2::FunctionSpaceTools fst; - MDFieldArrayFactory af("",ddims_,true); - - const panzer::PureBasis & basis = *(basis_layout->getBasis()); - - const int num_points = basis_layout->numPoints(); - const int num_basis = basis.cardinality(); - const int num_dim = basis_layout->dimension(); - const int num_cells = cub_points.extent(0); - - auto cell_basis_scalar = af.buildStaticArray("cell_basis_scalar",1,num_basis,num_points); - auto cell_cub_points = af.buildStaticArray("cell_cub_points",num_points,num_dim); - auto cell_grad_basis = af.buildStaticArray("cell_grad_basis",1,num_basis,num_points,num_dim); - auto cell_jac_inv = af.buildStaticArray("cell_jac_inv",1,num_points,num_dim,num_dim); - - auto cell_basis_ref_scalar = af.buildStaticArray("cell_basis_ref_scalar",num_basis,num_points); - auto cell_grad_basis_ref = af.buildStaticArray("cell_grad_basis_ref",num_basis,num_points,num_dim); - - for(int cell=0;cellgetValues(cell_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE); - - if(compute_derivatives){ - intrepid_basis->getValues(cell_grad_basis_ref.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_GRAD); - } - - // ============================================= - // Transform reference values to physical values - - fst::HGRADtransformVALUE(cell_basis_scalar.get_view(),cell_basis_ref_scalar.get_view()); - for(int b=0;b -void panzer::BasisValues2:: -evaluateValues_HCurl(const PHX::MDField & cub_points, - const PHX::MDField & jac, - const PHX::MDField & jac_det, - const PHX::MDField & jac_inv, - const PHX::MDField & weighted_measure) -{ - - TEUCHOS_ASSERT(getElementSpace() == PureBasis::HCURL); - - - typedef Intrepid2::FunctionSpaceTools fst; - MDFieldArrayFactory af("",ddims_,true); - - const panzer::PureBasis & basis = *(basis_layout->getBasis()); - - const int num_points = basis_layout->numPoints(); - const int num_basis = basis.cardinality(); - const int num_dim = basis_layout->dimension(); - const int num_cells = basis_layout->numCells(); - - auto cell_cub_points = af.buildStaticArray("cell_cub_points",num_points,num_dim); - auto cell_jac = af.buildStaticArray("cell_jac",1,num_points,num_dim,num_dim); - auto cell_jac_inv = af.buildStaticArray("cell_jac_inv",1,num_points,num_dim,num_dim); - auto cell_jac_det = af.buildStaticArray("cell_jac_det",1,num_points); - - auto cell_basis_vector = af.buildStaticArray("cell_basis_vector",1,num_basis,num_points,num_dim); - auto cell_curl_basis_scalar = af.buildStaticArray("cell_curl_basis_scalar",1,num_basis,num_points); - auto cell_curl_basis_vector = af.buildStaticArray("cell_curl_basis_vector",1,num_basis,num_points,num_dim); - - auto cell_curl_basis_ref = af.buildArray("cell_curl_basis_ref",num_basis,num_points,num_dim); - auto cell_curl_basis_ref_scalar = af.buildStaticArray("cell_curl_basis_ref_scalar",num_basis,num_points); - auto cell_basis_ref_vector = af.buildArray("cell_basis_ref_vector",num_basis,num_points,num_dim); - - for(int cell=0;cellgetValues(cell_basis_ref_vector.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE); - - if(compute_derivatives){ - if(num_dim==2){ - intrepid_basis->getValues(cell_curl_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_CURL); - } else if(num_dim==3){ - intrepid_basis->getValues(cell_curl_basis_ref.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_CURL); - } - } - - // ============================================= - // Transform reference values to physical values - - fst::HCURLtransformVALUE(cell_basis_vector.get_view(),cell_jac_inv.get_view(),cell_basis_ref_vector.get_view()); - for(int b=0;b -void panzer::BasisValues2:: -evaluateValues_HDiv(const PHX::MDField & cub_points, - const PHX::MDField & jac, - const PHX::MDField & jac_det, - const PHX::MDField & weighted_measure) -{ - - TEUCHOS_ASSERT(getElementSpace() == PureBasis::HDIV); - - typedef Intrepid2::FunctionSpaceTools fst; - MDFieldArrayFactory af("",ddims_,true); - - const panzer::PureBasis & basis = *(basis_layout->getBasis()); - - const int num_points = basis_layout->numPoints(); - const int num_basis = basis.cardinality(); - const int num_dim = basis_layout->dimension(); - const int num_cells = basis_layout->numCells(); - - auto cell_cub_points = af.buildStaticArray("cell_cub_points",num_points,num_dim); - auto cell_jac = af.buildStaticArray("cell_jac",1,num_points,num_dim,num_dim); - auto cell_jac_det = af.buildStaticArray("cell_jac_det",1,num_points); - - auto cell_basis_vector = af.buildStaticArray("cell_basis_vector",1,num_basis,num_points,num_dim); - auto cell_div_basis = af.buildStaticArray("cell_div_basis",1,num_basis,num_points); - - auto cell_basis_ref_vector = af.buildArray("cell_basis_ref_vector",num_basis,num_points,num_dim); - auto cell_div_basis_ref = af.buildStaticArray("cell_div_basis_ref",num_basis,num_points); - - for(int cell=0;cellgetValues(cell_basis_ref_vector.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE); - - if(compute_derivatives){ - intrepid_basis->getValues(cell_div_basis_ref.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_DIV); - } - - // ============================================= - // Transform reference values to physical values - - fst::HDIVtransformVALUE(cell_basis_vector.get_view(),cell_jac.get_view(),cell_jac_det.get_view(),cell_basis_ref_vector.get_view()); - for(int b=0;b -void panzer::BasisValues2:: -evaluateValuesCV(const PHX::MDField & cell_cub_points, - const PHX::MDField & jac, - const PHX::MDField & jac_det, - const PHX::MDField & jac_inv) -{ - MDFieldArrayFactory af("",ddims_,true); - - int num_ip = basis_layout->numPoints(); - int num_card = basis_layout->cardinality(); - int num_dim = basis_layout->dimension(); - - size_type num_cells = jac.extent(0); - - PureBasis::EElementSpace elmtspace = getElementSpace(); - ArrayDynamic dyn_cub_points = af.buildArray("dyn_cub_points", num_ip, num_dim); - - // Integration points are located on physical cells rather than reference cells, - // so we evaluate the basis in a loop over cells. - for (size_type icell = 0; icell < num_cells; ++icell) - { - for (int ip = 0; ip < num_ip; ++ip) - for (int d = 0; d < num_dim; ++d) - dyn_cub_points(ip,d) = cell_cub_points(icell,ip,d); - - if(elmtspace==PureBasis::CONST) { - ArrayDynamic dyn_basis_ref_scalar = af.buildArray("dyn_basis_ref_scalar",num_card,num_ip); - - intrepid_basis->getValues(dyn_basis_ref_scalar.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_VALUE); - - // transform values method just transfers values to array with cell index - no need to call - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_ip; ++ip) - basis_scalar(icell,b,ip) = dyn_basis_ref_scalar(b,ip); - - } - if(elmtspace==PureBasis::HGRAD) { - ArrayDynamic dyn_basis_ref_scalar = af.buildArray("dyn_basis_ref_scalar",num_card,num_ip); - - intrepid_basis->getValues(dyn_basis_ref_scalar.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_VALUE); - - // transform values method just transfers values to array with cell index - no need to call - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_ip; ++ip) - basis_scalar(icell,b,ip) = dyn_basis_ref_scalar(b,ip); - - if(compute_derivatives) { - - int one_cell = 1; - ArrayDynamic dyn_grad_basis_ref = af.buildArray("dyn_grad_basis_ref",num_card,num_ip,num_dim); - ArrayDynamic dyn_grad_basis = af.buildArray("dyn_grad_basis",one_cell,num_card,num_ip,num_dim); - ArrayDynamic dyn_jac_inv = af.buildArray("dyn_jac_inv",one_cell,num_ip,num_dim,num_dim); - - intrepid_basis->getValues(dyn_grad_basis_ref.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_GRAD); - - int cellInd = 0; - for (int ip = 0; ip < num_ip; ++ip) - for (int d1 = 0; d1 < num_dim; ++d1) - for (int d2 = 0; d2 < num_dim; ++d2) - dyn_jac_inv(cellInd,ip,d1,d2) = jac_inv(icell,ip,d1,d2); - - Intrepid2::FunctionSpaceTools::HGRADtransformGRAD(dyn_grad_basis.get_view(), - dyn_jac_inv.get_view(), - dyn_grad_basis_ref.get_view()); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_ip; ++ip) - for (int d = 0; d < num_dim; ++d) - grad_basis(icell,b,ip,d) = dyn_grad_basis(0,b,ip,d); - - } - } - else if(elmtspace==PureBasis::HCURL) { - ArrayDynamic dyn_basis_ref_vector = af.buildArray("dyn_basis_ref_vector",num_card,num_ip,num_dim); - - intrepid_basis->getValues(dyn_basis_ref_vector.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_VALUE); - - int one_cell = 1; - ArrayDynamic dyn_basis_vector = af.buildArray("dyn_basis_vector",one_cell,num_card,num_ip,num_dim); - ArrayDynamic dyn_jac_inv = af.buildArray("dyn_jac_inv",one_cell,num_ip,num_dim,num_dim); - - int cellInd = 0; - for (int ip = 0; ip < num_ip; ++ip) - for (int d1 = 0; d1 < num_dim; ++d1) - for (int d2 = 0; d2 < num_dim; ++d2) - dyn_jac_inv(cellInd,ip,d1,d2) = jac_inv(icell,ip,d1,d2); - - Intrepid2::FunctionSpaceTools::HCURLtransformVALUE(dyn_basis_vector.get_view(), - dyn_jac_inv.get_view(), - dyn_basis_ref_vector.get_view()); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_ip; ++ip) - for (int d = 0; d < num_dim; ++d) - basis_vector(icell,b,ip,d) = dyn_basis_vector(0,b,ip,d); - - if(compute_derivatives && num_dim ==2) { - - int one_cell = 1; - ArrayDynamic dyn_curl_basis_ref_scalar = af.buildArray("dyn_curl_basis_ref_scalar",num_card,num_ip); - ArrayDynamic dyn_curl_basis_scalar = af.buildArray("dyn_curl_basis_scalar",one_cell,num_card,num_ip); - ArrayDynamic dyn_jac_det = af.buildArray("dyn_jac_det",one_cell,num_ip); - - intrepid_basis->getValues(dyn_curl_basis_ref_scalar.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_CURL); - - int cellInd = 0; - for (int ip = 0; ip < num_ip; ++ip) - dyn_jac_det(cellInd,ip) = jac_det(icell,ip); - - Intrepid2::FunctionSpaceTools::HDIVtransformDIV(dyn_curl_basis_scalar.get_view(), - dyn_jac_det.get_view(), - dyn_curl_basis_ref_scalar.get_view()); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_ip; ++ip) - curl_basis_scalar(icell,b,ip) = dyn_curl_basis_scalar(0,b,ip); - - } - if(compute_derivatives && num_dim ==3) { - - int one_cell = 1; - ArrayDynamic dyn_curl_basis_ref = af.buildArray("dyn_curl_basis_ref_vector",num_card,num_ip,num_dim); - ArrayDynamic dyn_curl_basis = af.buildArray("dyn_curl_basis_vector",one_cell,num_card,num_ip,num_dim); - ArrayDynamic dyn_jac_det = af.buildArray("dyn_jac_det",one_cell,num_ip); - ArrayDynamic dyn_jac = af.buildArray("dyn_jac",one_cell,num_ip,num_dim,num_dim); - - intrepid_basis->getValues(dyn_curl_basis_ref.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_CURL); - - int cellInd = 0; - for (int ip = 0; ip < num_ip; ++ip) - { - dyn_jac_det(cellInd,ip) = jac_det(icell,ip); - for (int d1 = 0; d1 < num_dim; ++d1) - for (int d2 = 0; d2 < num_dim; ++d2) - dyn_jac(cellInd,ip,d1,d2) = jac(icell,ip,d1,d2); - } - - Intrepid2::FunctionSpaceTools::HCURLtransformCURL(dyn_curl_basis.get_view(), - dyn_jac.get_view(), - dyn_jac_det.get_view(), - dyn_curl_basis_ref.get_view()); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_ip; ++ip) - for (int d = 0; d < num_dim; ++d) - curl_basis_vector(icell,b,ip,d) = dyn_curl_basis(0,b,ip,d); - - } - - } - else if(elmtspace==PureBasis::HDIV) { - - ArrayDynamic dyn_basis_ref_vector = af.buildArray("dyn_basis_ref_vector",num_card,num_ip,num_dim); - - intrepid_basis->getValues(dyn_basis_ref_vector.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_VALUE); - - int one_cell= 1; - ArrayDynamic dyn_basis_vector = af.buildArray("dyn_basis_vector",one_cell,num_card,num_ip,num_dim); - ArrayDynamic dyn_jac = af.buildArray("dyn_jac",one_cell,num_ip,num_dim,num_dim); - ArrayDynamic dyn_jac_det = af.buildArray("dyn_jac_det",one_cell,num_ip); - - int cellInd = 0; - for (int ip = 0; ip < num_ip; ++ip) - { - dyn_jac_det(cellInd,ip) = jac_det(icell,ip); - for (int d1 = 0; d1 < num_dim; ++d1) - for (int d2 = 0; d2 < num_dim; ++d2) - dyn_jac(cellInd,ip,d1,d2) = jac(icell,ip,d1,d2); - } - - Intrepid2::FunctionSpaceTools::HDIVtransformVALUE(dyn_basis_vector.get_view(), - dyn_jac.get_view(), - dyn_jac_det.get_view(), - dyn_basis_ref_vector.get_view()); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_ip; ++ip) - for (int d = 0; d < num_dim; ++d) - basis_vector(icell,b,ip,d) = dyn_basis_vector(0,b,ip,d); - - if(compute_derivatives) { - - ArrayDynamic dyn_div_basis_ref = af.buildArray("dyn_div_basis_ref_scalar",num_card,num_ip); - ArrayDynamic dyn_div_basis = af.buildArray("dyn_div_basis_scalar",one_cell,num_card,num_ip); - - intrepid_basis->getValues(dyn_div_basis_ref.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_DIV); - - Intrepid2::FunctionSpaceTools::HDIVtransformDIV(dyn_div_basis.get_view(), - dyn_jac_det.get_view(), - dyn_div_basis_ref.get_view()); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_ip; ++ip) - div_basis(icell,b,ip) = dyn_div_basis(0,b,ip); - - } - - } - else { TEUCHOS_ASSERT(false); } - - } // cell loop - -} - -template -void panzer::BasisValues2:: -evaluateReferenceValues(const PHX::MDField & cub_points,bool compute_derivatives,bool use_vertex_coordinates) -{ - MDFieldArrayFactory af("",ddims_,true); - - int num_quad = basis_layout->numPoints(); - int num_dim = basis_layout->dimension(); - int num_card = basis_layout->cardinality(); - - ArrayDynamic dyn_cub_points = af.buildArray("dyn_cub_points", num_quad,num_dim); - - for (int ip = 0; ip < num_quad; ++ip) - for (int d = 0; d < num_dim; ++d) - dyn_cub_points(ip,d) = cub_points(ip,d); - - PureBasis::EElementSpace elmtspace = getElementSpace(); - if(elmtspace==PureBasis::HGRAD || elmtspace==PureBasis::CONST) { - ArrayDynamic dyn_basis_ref_scalar = af.buildArray("dyn_basis_ref_scalar",num_card,num_quad); - - intrepid_basis->getValues(dyn_basis_ref_scalar.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_VALUE); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_quad; ++ip) - basis_ref_scalar(b,ip) = dyn_basis_ref_scalar(b,ip); - } - else if(elmtspace==PureBasis::HDIV || elmtspace==PureBasis::HCURL) { - ArrayDynamic dyn_basis_ref_vector = af.buildArray("dyn_basis_ref_vector",num_card,num_quad,num_dim); - - intrepid_basis->getValues(dyn_basis_ref_vector.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_VALUE); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_quad; ++ip) - for (int d = 0; d < num_dim; ++d) - basis_ref_vector(b,ip,d) = dyn_basis_ref_vector(b,ip,d); - } - else { TEUCHOS_ASSERT(false); } - - if(elmtspace==PureBasis::HGRAD && compute_derivatives) { - ArrayDynamic dyn_grad_basis_ref = af.buildArray("dyn_basis_ref_vector",num_card,num_quad,num_dim); - - intrepid_basis->getValues(dyn_grad_basis_ref.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_GRAD); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_quad; ++ip) - for (int d = 0; d < num_dim; ++d) - grad_basis_ref(b,ip,d) = dyn_grad_basis_ref(b,ip,d); - } - else if(elmtspace==PureBasis::HCURL && compute_derivatives && num_dim==2) { - ArrayDynamic dyn_curl_basis_ref = af.buildArray("dyn_curl_basis_ref_scalar",num_card,num_quad); - - intrepid_basis->getValues(dyn_curl_basis_ref.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_CURL); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_quad; ++ip) - curl_basis_ref_scalar(b,ip) = dyn_curl_basis_ref(b,ip); - } - else if(elmtspace==PureBasis::HCURL && compute_derivatives && num_dim==3) { - ArrayDynamic dyn_curl_basis_ref = af.buildArray("dyn_curl_basis_ref_vector",num_card,num_quad,num_dim); - - intrepid_basis->getValues(dyn_curl_basis_ref.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_CURL); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_quad; ++ip) - for (int d = 0; d < num_dim; ++d) - curl_basis_ref_vector(b,ip,d) = dyn_curl_basis_ref(b,ip,d); - } - else if(elmtspace==PureBasis::HDIV && compute_derivatives) { - ArrayDynamic dyn_div_basis_ref = af.buildArray("dyn_div_basis_ref_scalar",num_card,num_quad); - - intrepid_basis->getValues(dyn_div_basis_ref.get_view(), - dyn_cub_points.get_view(), - Intrepid2::OPERATOR_DIV); - - for (int b = 0; b < num_card; ++b) - for (int ip = 0; ip < num_quad; ++ip) - div_basis_ref(b,ip) = dyn_div_basis_ref(b,ip); - } - - - if(use_vertex_coordinates) { - // Intrepid removes fad types from the coordinate scalar type. We - // pull the actual field scalar type from the basis object to be - // consistent. - if (elmtspace != PureBasis::CONST) { - using coordsScalarType = typename Intrepid2::Basis::scalarType; - auto dyn_basis_coordinates_ref = af.buildArray("basis_coordinates_ref", - basis_coordinates_ref.extent(0), - basis_coordinates_ref.extent(1)); - intrepid_basis->getDofCoords(dyn_basis_coordinates_ref.get_view()); - - // fill in basis coordinates - for (int i = 0; i < basis_coordinates_ref.extent_int(0); ++i) - for (int j = 0; j < basis_coordinates_ref.extent_int(1); ++j) - basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j); - } - } - - references_evaluated = true; -} - -// method for applying orientations -template -void BasisValues2:: -applyOrientations(const std::vector & orientations) -{ - if (!intrepid_basis->requireOrientation()) - return; - - typedef Intrepid2::OrientationTools ots; - const PureBasis::EElementSpace elmtspace = getElementSpace(); - - // orientation (right now std vector) is created using push_back method. - // thus, its size is the actual number of elements to be applied. - // on the other hand, basis_layout num cell indicates workset size. - // to get right size of cells, use minimum of them. - const int num_cell_basis_layout = basis_layout->numCells(), num_cell_orientation = orientations.size(); - const int num_cell = num_cell_basis_layout < num_cell_orientation ? num_cell_basis_layout : num_cell_orientation; - const int num_dim = basis_layout->dimension(); - const Kokkos::pair range_cell(0, num_cell); - - Kokkos::DynRankView - drv_orts((Intrepid2::Orientation*)orientations.data(), num_cell); - - /// - /// HGRAD elements - /// - if (elmtspace==PureBasis::HGRAD) { - { - { - auto drv_basis_scalar = Kokkos::subview(basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); - auto drv_basis_scalar_tmp = Kokkos::createDynRankView(basis_scalar.get_view(), - "drv_basis_scalar_tmp", - drv_basis_scalar.extent(0), // C - drv_basis_scalar.extent(1), // F - drv_basis_scalar.extent(2)); // P - Kokkos::deep_copy(drv_basis_scalar_tmp, drv_basis_scalar); - ots::modifyBasisByOrientation(drv_basis_scalar, - drv_basis_scalar_tmp, - drv_orts, - intrepid_basis); - } - if(build_weighted) { - auto drv_basis_scalar = Kokkos::subview(weighted_basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); - auto drv_basis_scalar_tmp = Kokkos::createDynRankView(weighted_basis_scalar.get_view(), - "drv_basis_scalar_tmp", - drv_basis_scalar.extent(0), // C - drv_basis_scalar.extent(1), // F - drv_basis_scalar.extent(2)); // P - Kokkos::deep_copy(drv_basis_scalar_tmp, drv_basis_scalar); - ots::modifyBasisByOrientation(drv_basis_scalar, - drv_basis_scalar_tmp, - drv_orts, - intrepid_basis); - } - - } - - if (compute_derivatives) { - { - auto drv_grad_basis = Kokkos::subview(grad_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - auto drv_grad_basis_tmp = Kokkos::createDynRankView(grad_basis.get_view(), - "drv_grad_basis_tmp", - drv_grad_basis.extent(0), // C - drv_grad_basis.extent(1), // F - drv_grad_basis.extent(2), // P - drv_grad_basis.extent(3)); // D - Kokkos::deep_copy(drv_grad_basis_tmp, drv_grad_basis); - ots::modifyBasisByOrientation(drv_grad_basis, - drv_grad_basis_tmp, - drv_orts, - intrepid_basis); - } - if(build_weighted) { - auto drv_grad_basis = Kokkos::subview(weighted_grad_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - auto drv_grad_basis_tmp = Kokkos::createDynRankView(weighted_grad_basis.get_view(), - "drv_grad_basis_tmp", - drv_grad_basis.extent(0), // C - drv_grad_basis.extent(1), // F - drv_grad_basis.extent(2), // P - drv_grad_basis.extent(3)); // D - Kokkos::deep_copy(drv_grad_basis_tmp, drv_grad_basis); - ots::modifyBasisByOrientation(drv_grad_basis, - drv_grad_basis_tmp, - drv_orts, - intrepid_basis); - } - } - } - - /// - /// hcurl 2d elements - /// - else if (elmtspace==PureBasis::HCURL && num_dim==2) { - { - { - auto drv_basis_vector = Kokkos::subview(basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - auto drv_basis_vector_tmp = Kokkos::createDynRankView(basis_vector.get_view(), - "drv_basis_vector_tmp", - drv_basis_vector.extent(0), // C - drv_basis_vector.extent(1), // F - drv_basis_vector.extent(2), // P - drv_basis_vector.extent(3)); // D - Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); - ots::modifyBasisByOrientation(drv_basis_vector, - drv_basis_vector_tmp, - drv_orts, - intrepid_basis); - } - if(build_weighted) { - auto drv_basis_vector = Kokkos::subview(weighted_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - auto drv_basis_vector_tmp = Kokkos::createDynRankView(weighted_basis_vector.get_view(), - "drv_basis_vector_tmp", - drv_basis_vector.extent(0), // C - drv_basis_vector.extent(1), // F - drv_basis_vector.extent(2), // P - drv_basis_vector.extent(3)); // D - Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); - ots::modifyBasisByOrientation(drv_basis_vector, - drv_basis_vector_tmp, - drv_orts, - intrepid_basis); - } - } - - if (compute_derivatives) { - { - auto drv_curl_basis_scalar = Kokkos::subview(curl_basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); - auto drv_curl_basis_scalar_tmp = Kokkos::createDynRankView(curl_basis_scalar.get_view(), - "drv_curl_basis_scalar_tmp", - drv_curl_basis_scalar.extent(0), // C - drv_curl_basis_scalar.extent(1), // F - drv_curl_basis_scalar.extent(2)); // P - Kokkos::deep_copy(drv_curl_basis_scalar_tmp, drv_curl_basis_scalar); - ots::modifyBasisByOrientation(drv_curl_basis_scalar, - drv_curl_basis_scalar_tmp, - drv_orts, - intrepid_basis); - } - - if(build_weighted) { - auto drv_curl_basis_scalar = Kokkos::subview(weighted_curl_basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); - auto drv_curl_basis_scalar_tmp = Kokkos::createDynRankView(weighted_curl_basis_scalar.get_view(), - "drv_curl_basis_scalar_tmp", - drv_curl_basis_scalar.extent(0), // C - drv_curl_basis_scalar.extent(1), // F - drv_curl_basis_scalar.extent(2)); // P - Kokkos::deep_copy(drv_curl_basis_scalar_tmp, drv_curl_basis_scalar); - ots::modifyBasisByOrientation(drv_curl_basis_scalar, - drv_curl_basis_scalar_tmp, - drv_orts, - intrepid_basis); - } - } - } - - /// - /// hcurl 3d elements - /// - else if (elmtspace==PureBasis::HCURL && num_dim==3) { - { - { - auto drv_basis_vector = Kokkos::subview(basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - auto drv_basis_vector_tmp = Kokkos::createDynRankView(basis_vector.get_view(), - "drv_basis_vector_tmp", - drv_basis_vector.extent(0), // C - drv_basis_vector.extent(1), // F - drv_basis_vector.extent(2), // P - drv_basis_vector.extent(3)); // D - Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); - ots::modifyBasisByOrientation(drv_basis_vector, - drv_basis_vector_tmp, - drv_orts, - intrepid_basis); - } - if(build_weighted) { - auto drv_basis_vector = Kokkos::subview(weighted_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - auto drv_basis_vector_tmp = Kokkos::createDynRankView(weighted_basis_vector.get_view(), - "drv_basis_vector_tmp", - drv_basis_vector.extent(0), // C - drv_basis_vector.extent(1), // F - drv_basis_vector.extent(2), // P - drv_basis_vector.extent(3)); // D - Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); - ots::modifyBasisByOrientation(drv_basis_vector, - drv_basis_vector_tmp, - drv_orts, - intrepid_basis); - } - } - - if (compute_derivatives) { - { - auto drv_curl_basis_vector = Kokkos::subview(curl_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - auto drv_curl_basis_vector_tmp = Kokkos::createDynRankView(curl_basis_vector.get_view(), - "drv_curl_basis_vector_tmp", - drv_curl_basis_vector.extent(0), // C - drv_curl_basis_vector.extent(1), // F - drv_curl_basis_vector.extent(2), // P - drv_curl_basis_vector.extent(3)); // D - Kokkos::deep_copy(drv_curl_basis_vector_tmp, drv_curl_basis_vector); - ots::modifyBasisByOrientation(drv_curl_basis_vector, - drv_curl_basis_vector_tmp, - drv_orts, - intrepid_basis); - } - if(build_weighted) { - auto drv_curl_basis_vector = Kokkos::subview(weighted_curl_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - auto drv_curl_basis_vector_tmp = Kokkos::createDynRankView(weighted_curl_basis_vector.get_view(), - "drv_curl_basis_vector_tmp", - drv_curl_basis_vector.extent(0), // C - drv_curl_basis_vector.extent(1), // F - drv_curl_basis_vector.extent(2), // P - drv_curl_basis_vector.extent(3)); // D - Kokkos::deep_copy(drv_curl_basis_vector_tmp, drv_curl_basis_vector); - ots::modifyBasisByOrientation(drv_curl_basis_vector, - drv_curl_basis_vector_tmp, - drv_orts, - intrepid_basis); - } - } - } - /// - /// hdiv elements (2d and 3d) - /// - else if (elmtspace==PureBasis::HDIV) { - { - { - auto drv_basis_vector = Kokkos::subview(basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - auto drv_basis_vector_tmp = Kokkos::createDynRankView(basis_vector.get_view(), - "drv_basis_vector_tmp", - drv_basis_vector.extent(0), // C - drv_basis_vector.extent(1), // F - drv_basis_vector.extent(2), // P - drv_basis_vector.extent(3)); // D - Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); - ots::modifyBasisByOrientation(drv_basis_vector, - drv_basis_vector_tmp, - drv_orts, - intrepid_basis); - } - if(build_weighted) { - auto drv_basis_vector = Kokkos::subview(weighted_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - auto drv_basis_vector_tmp = Kokkos::createDynRankView(weighted_basis_vector.get_view(), - "drv_basis_vector_tmp", - drv_basis_vector.extent(0), // C - drv_basis_vector.extent(1), // F - drv_basis_vector.extent(2), // P - drv_basis_vector.extent(3)); // D - Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); - ots::modifyBasisByOrientation(drv_basis_vector, - drv_basis_vector_tmp, - drv_orts, - intrepid_basis); - } - } - if (compute_derivatives) { - { - auto drv_div_basis = Kokkos::subview(div_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); - auto drv_div_basis_tmp = Kokkos::createDynRankView(div_basis.get_view(), - "drv_div_basis_tmp", - drv_div_basis.extent(0), // C - drv_div_basis.extent(1), // F - drv_div_basis.extent(2)); // P - Kokkos::deep_copy(drv_div_basis_tmp, drv_div_basis); - ots::modifyBasisByOrientation(drv_div_basis, - drv_div_basis_tmp, - drv_orts, - intrepid_basis); - } - if(build_weighted) { - auto drv_div_basis = Kokkos::subview(weighted_div_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); - auto drv_div_basis_tmp = Kokkos::createDynRankView(weighted_div_basis.get_view(), - "drv_div_basis_tmp", - drv_div_basis.extent(0), // C - drv_div_basis.extent(1), // F - drv_div_basis.extent(2)); // P - Kokkos::deep_copy(drv_div_basis_tmp, drv_div_basis); - ots::modifyBasisByOrientation(drv_div_basis, - drv_div_basis_tmp, - drv_orts, - intrepid_basis); - } - } - } -} - -// method for applying orientations -template -void BasisValues2:: -applyOrientations(const PHX::MDField & orientations) -{ - - TEUCHOS_TEST_FOR_EXCEPT_MSG(true,"panzer::BasisValues2::applyOrientations : this should not be called."); - - int num_cell = orientations.extent(0); - int num_basis = orientations.extent(1); - int num_dim = basis_layout->dimension(); - int num_ip = basis_layout->numPoints(); - PureBasis::EElementSpace elmtspace = getElementSpace(); - - if(elmtspace==PureBasis::HCURL && num_dim==2) { - - // setup the orientations for the trial space - // Intrepid2::FunctionSpaceTools::applyFieldSigns(basis_vector,orientations); - - for (int c=0; c(curl_basis_scalar,orientations); - for (int c=0; c::applyFieldSigns(weighted_basis_vector.get_view(),orientations.get_view()); - - if(compute_derivatives) - Intrepid2::FunctionSpaceTools::applyFieldSigns(weighted_curl_basis_scalar.get_view(),orientations.get_view()); - } - } - else if(elmtspace==PureBasis::HCURL && num_dim==3) { - - // setup the orientations for the trial space - // Intrepid2::FunctionSpaceTools::applyFieldSigns(basis_vector,orientations); - - for (int c=0; c(curl_basis_vector,orientations); - for (int c=0; c::applyFieldSigns(weighted_basis_vector.get_view(),orientations.get_view()); - - if(compute_derivatives) - Intrepid2::FunctionSpaceTools::applyFieldSigns(weighted_curl_basis_vector.get_view(),orientations.get_view()); - } - } - else if(elmtspace==PureBasis::HDIV) { - // setup the orientations for the trial space - // Intrepid2::FunctionSpaceTools::applyFieldSigns(basis_vector,orientations); - - for (int c=0; c(div_basis,orientations); - - for (int c=0; c::applyFieldSigns(weighted_basis_vector.get_view(),orientations.get_view()); - - if(compute_derivatives) - Intrepid2::FunctionSpaceTools::applyFieldSigns(weighted_div_basis.get_view(),orientations.get_view()); - } - } -} - -template -PureBasis::EElementSpace BasisValues2::getElementSpace() const -{ return basis_layout->getBasis()->getElementSpace(); } - -template -void panzer::BasisValues2:: -setupArrays(const Teuchos::RCP& layout, - bool computeDerivatives) -{ - MDFieldArrayFactory af(prefix,alloc_arrays); - - compute_derivatives = computeDerivatives; - basis_layout = layout; - Teuchos::RCP basisDesc = layout->getBasis(); - - // for convience pull out basis and quadrature information - int num_quad = layout->numPoints(); - int dim = basisDesc->dimension(); - int card = basisDesc->cardinality(); - int numcells = basisDesc->numCells(); - panzer::PureBasis::EElementSpace elmtspace = basisDesc->getElementSpace(); - Teuchos::RCP cellTopo = basisDesc->getCellTopology(); - - intrepid_basis = basisDesc->getIntrepid2Basis(); - - // allocate field containers - // field sizes defined by http://trilinos.sandia.gov/packages/docs/dev/packages/intrepid/doc/html/basis_page.html#basis_md_array_sec - - // compute basis fields - if(elmtspace==panzer::PureBasis::HGRAD) { - // HGRAD is a nodal field - - // build values - /////////////////////////////////////////////// - basis_ref_scalar = af.buildStaticArray("basis_ref",card,num_quad); // F, P - basis_scalar = af.buildStaticArray("basis",numcells,card,num_quad); - - if(build_weighted) - weighted_basis_scalar = af.buildStaticArray("weighted_basis",numcells,card,num_quad); - - // build gradients - /////////////////////////////////////////////// - - if(compute_derivatives) { - grad_basis_ref = af.buildStaticArray("grad_basis_ref",card,num_quad,dim); // F, P, D - grad_basis = af.buildStaticArray("grad_basis",numcells,card,num_quad,dim); - - if(build_weighted) - weighted_grad_basis = af.buildStaticArray("weighted_grad_basis",numcells,card,num_quad,dim); - } - - // build curl - /////////////////////////////////////////////// - - // nothing - HGRAD does not support CURL operation - } - else if(elmtspace==panzer::PureBasis::HCURL) { - // HCURL is a vector field - - // build values - /////////////////////////////////////////////// - - basis_ref_vector = af.buildStaticArray("basis_ref",card,num_quad,dim); // F, P, D - basis_vector = af.buildStaticArray("basis",numcells,card,num_quad,dim); - - if(build_weighted) - weighted_basis_vector = af.buildStaticArray("weighted_basis",numcells,card,num_quad,dim); - - // build gradients - /////////////////////////////////////////////// - - // nothing - HCURL does not support GRAD operation - - // build curl - /////////////////////////////////////////////// - - if(compute_derivatives) { - if(dim==2) { - // curl of HCURL basis is not dimension dependent - curl_basis_ref_scalar = af.buildStaticArray("curl_basis_ref",card,num_quad); // F, P - curl_basis_scalar = af.buildStaticArray("curl_basis",numcells,card,num_quad); - - if(build_weighted) - weighted_curl_basis_scalar = af.buildStaticArray("weighted_curl_basis",numcells,card,num_quad); - } - else if(dim==3){ - curl_basis_ref_vector = af.buildStaticArray("curl_basis_ref",card,num_quad,dim); // F, P, D - curl_basis_vector = af.buildStaticArray("curl_basis",numcells,card,num_quad,dim); - - if(build_weighted) - weighted_curl_basis_vector = af.buildStaticArray("weighted_curl_basis",numcells,card,num_quad,dim); - } - else { TEUCHOS_ASSERT(false); } // what do I do with 1D? - } - } - else if(elmtspace==panzer::PureBasis::HDIV) { - // HDIV is a vector field - - // build values - /////////////////////////////////////////////// - - basis_ref_vector = af.buildStaticArray("basis_ref",card,num_quad,dim); // F, P, D - basis_vector = af.buildStaticArray("basis",numcells,card,num_quad,dim); - - if(build_weighted) - weighted_basis_vector = af.buildStaticArray("weighted_basis",numcells,card,num_quad,dim); - - // build gradients - /////////////////////////////////////////////// - - // nothing - HCURL does not support GRAD operation - - // build curl - /////////////////////////////////////////////// - - // nothing - HDIV does not support CURL operation - - // build div - /////////////////////////////////////////////// - - if(compute_derivatives) { - div_basis_ref = af.buildStaticArray("div_basis_ref",card,num_quad); // F, P - div_basis = af.buildStaticArray("div_basis",numcells,card,num_quad); - - if(build_weighted) - weighted_div_basis = af.buildStaticArray("weighted_div_basis",numcells,card,num_quad); - } - } - else if(elmtspace==panzer::PureBasis::CONST) { - // CONST is a nodal field - - // build values - /////////////////////////////////////////////// - basis_ref_scalar = af.buildStaticArray("basis_ref",card,num_quad); // F, P - basis_scalar = af.buildStaticArray("basis",numcells,card,num_quad); - - if(build_weighted) - weighted_basis_scalar = af.buildStaticArray("weighted_basis",numcells,card,num_quad); - - // build gradients - /////////////////////////////////////////////// - - // nothing - CONST does not support GRAD operation - - // build curl - /////////////////////////////////////////////// - - // nothing - CONST does not support CURL operation - - // build div - /////////////////////////////////////////////// - - // nothing - CONST does not support DIV operation - } - else { TEUCHOS_ASSERT(false); } - - basis_coordinates_ref = af.buildStaticArray("basis_coordinates_ref",card,dim); - basis_coordinates = af.buildStaticArray("basis_coordinates",numcells,card,dim); -} - -// do some explicit instantiation so things compile faster. - #define BASIS_VALUES_INSTANTIATION(SCALAR) \ template class BasisValues2; BASIS_VALUES_INSTANTIATION(panzer::Traits::RealType) -BASIS_VALUES_INSTANTIATION(panzer::Traits::FadType) + +// Disabled due to long build times on cuda (30+ minutes for this +// instantiaiton alone. If we need sensitivities wrt coordinates, we +// can reenable. +//BASIS_VALUES_INSTANTIATION(panzer::Traits::FadType) + #ifdef Panzer_BUILD_HESSIAN_SUPPORT BASIS_VALUES_INSTANTIATION(panzer::Traits::HessianType) #endif -} // namespace panzer +} diff --git a/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp b/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp new file mode 100644 index 000000000000..f07ddfe6cfd5 --- /dev/null +++ b/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp @@ -0,0 +1,1548 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// 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 Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#include "PanzerDiscFE_config.hpp" +#include "Panzer_Traits.hpp" + +#include "Panzer_CommonArrayFactories.hpp" +#include "Kokkos_ViewFactory.hpp" + +#include "Intrepid2_Utils.hpp" +#include "Intrepid2_FunctionSpaceTools.hpp" +#include "Intrepid2_Orientation.hpp" +#include "Intrepid2_OrientationTools.hpp" + + +namespace panzer { + +template +void panzer::BasisValues2:: +evaluateValues(const PHX::MDField & cub_points, + const PHX::MDField & jac, + const PHX::MDField & jac_det, + const PHX::MDField & jac_inv) +{ + PHX::MDField weighted_measure; + PHX::MDField vertex_coordinates; + build_weighted = false; + evaluateValues(cub_points,jac,jac_det,jac_inv,weighted_measure,vertex_coordinates,false); +} + +template +void panzer::BasisValues2:: +evaluateValues(const PHX::MDField & cub_points, + const PHX::MDField & jac, + const PHX::MDField & jac_det, + const PHX::MDField & jac_inv, + const PHX::MDField & weighted_measure, + const PHX::MDField & vertex_coordinates, + bool use_vertex_coordinates) +{ + MDFieldArrayFactory af("",ddims_,true); + + int num_dim = basis_layout->dimension(); + + // currently this just copies on the basis objects are converted + // in intrepid + evaluateReferenceValues(cub_points,compute_derivatives,use_vertex_coordinates); + + PureBasis::EElementSpace elmtspace = getElementSpace(); + if(elmtspace==PureBasis::CONST || + elmtspace==PureBasis::HGRAD) { + Intrepid2::FunctionSpaceTools:: + HGRADtransformVALUE(basis_scalar.get_view(), + basis_ref_scalar.get_view()); + + if(build_weighted) { + Intrepid2::FunctionSpaceTools:: + multiplyMeasure(weighted_basis_scalar.get_view(), + weighted_measure.get_view(), + basis_scalar.get_view()); + } + } + else if(elmtspace==PureBasis::HCURL) { + Intrepid2::FunctionSpaceTools:: + HCURLtransformVALUE(basis_vector.get_view(), + jac_inv.get_view(), + basis_ref_vector.get_view()); + + if(build_weighted) { + Intrepid2::FunctionSpaceTools:: + multiplyMeasure(weighted_basis_vector.get_view(), + weighted_measure.get_view(), + basis_vector.get_view()); + } + } + else if(elmtspace==PureBasis::HDIV) + { + Intrepid2::FunctionSpaceTools:: + HDIVtransformVALUE(basis_vector.get_view(), + jac.get_view(), + jac_det.get_view(), + basis_ref_vector.get_view()); + + if(build_weighted) { + Intrepid2::FunctionSpaceTools:: + multiplyMeasure(weighted_basis_vector.get_view(), + weighted_measure.get_view(), + basis_vector.get_view()); + } + } + else { TEUCHOS_ASSERT(false); } + + if(elmtspace==PureBasis::HGRAD && compute_derivatives) { + Intrepid2::FunctionSpaceTools:: + HGRADtransformGRAD(grad_basis.get_view(), + jac_inv.get_view(), + grad_basis_ref.get_view()); + + if(build_weighted) { + Intrepid2::FunctionSpaceTools:: + multiplyMeasure(weighted_grad_basis.get_view(), + weighted_measure.get_view(), + grad_basis.get_view()); + } + } + else if(elmtspace==PureBasis::HCURL && num_dim==2 && compute_derivatives) { + Intrepid2::FunctionSpaceTools:: + HDIVtransformDIV(curl_basis_scalar.get_view(), + jac_det.get_view(), // note only volume deformation is needed! + // this relates directly to this being in + // the divergence space in 2D! + curl_basis_ref_scalar.get_view()); + + if(build_weighted) { + Intrepid2::FunctionSpaceTools:: + multiplyMeasure(weighted_curl_basis_scalar.get_view(), + weighted_measure.get_view(), + curl_basis_scalar.get_view()); + } + } + else if(elmtspace==PureBasis::HCURL && num_dim==3 && compute_derivatives) { + Intrepid2::FunctionSpaceTools:: + HCURLtransformCURL(curl_basis_vector.get_view(), + jac.get_view(), + jac_det.get_view(), + curl_basis_ref_vector.get_view()); + + if(build_weighted) { + Intrepid2::FunctionSpaceTools:: + multiplyMeasure(weighted_curl_basis_vector.get_view(), + weighted_measure.get_view(), + curl_basis_vector.get_view()); + } + } + else if(elmtspace==PureBasis::HDIV && compute_derivatives) { + Intrepid2::FunctionSpaceTools:: + HDIVtransformDIV(div_basis.get_view(), + jac_det.get_view(), + div_basis_ref.get_view()); + + if(build_weighted) { + Intrepid2::FunctionSpaceTools:: + multiplyMeasure(weighted_div_basis.get_view(), + weighted_measure.get_view(), + div_basis.get_view()); + } + } + + // If basis supports coordinate values at basis points, then + // compute these values + if(use_vertex_coordinates) { + // Teuchos::RCP > coords + // = Teuchos::rcp_dynamic_cast >(intrepid_basis); + // if (!Teuchos::is_null(coords)) { +/* + ArrayDynamic dyn_basis_coordinates_ref = af.buildArray("basis_coordinates_ref",basis_coordinates_ref.extent(0),basis_coordinates_ref.extent(1)); + coords->getDofCoords(dyn_basis_coordinates_ref); + + // fill in basis coordinates + for (size_type i = 0; i < basis_coordinates_ref.extent(0); ++i) + for (size_type j = 0; j < basis_coordinates_ref.extent(1); ++j) + basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j); +*/ + + Intrepid2::CellTools cell_tools; + cell_tools.mapToPhysicalFrame(basis_coordinates.get_view(), + basis_coordinates_ref.get_view(), + vertex_coordinates.get_view(), + intrepid_basis->getBaseCellTopology()); + } +} + + + + + + + +template +void panzer::BasisValues2:: +evaluateBasisCoordinates(const PHX::MDField & vertex_coordinates) +{ + MDFieldArrayFactory af("",ddims_,true); + + // intrepid_basis->getDofCoords requires DynRankView, but basis_coordinates_ref is more of a static View + // We use an auxiliary 'dyn' array to get around this + using coordsScalarType = typename Intrepid2::Basis::scalarType; + auto dyn_basis_coordinates_ref = af.buildArray("basis_coordinates_ref", + basis_coordinates_ref.extent(0), + basis_coordinates_ref.extent(1)); + intrepid_basis->getDofCoords(dyn_basis_coordinates_ref.get_view()); + + // fill in basis coordinates + for (int i = 0; i < basis_coordinates_ref.extent_int(0); ++i) + for (int j = 0; j < basis_coordinates_ref.extent_int(1); ++j) + basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j); + + Intrepid2::CellTools cell_tools; + cell_tools.mapToPhysicalFrame(basis_coordinates.get_view(), + basis_coordinates_ref.get_view(), + vertex_coordinates.get_view(), + intrepid_basis->getBaseCellTopology()); +} + + + + + +template +void panzer::BasisValues2:: +evaluateValues(const PHX::MDField & cub_points, + const PHX::MDField & jac, + const PHX::MDField & jac_det, + const PHX::MDField & jac_inv, + const PHX::MDField & weighted_measure, + const PHX::MDField & vertex_coordinates, + bool use_vertex_coordinates) +{ + + PureBasis::EElementSpace elmtspace = getElementSpace(); + + if(elmtspace == PureBasis::CONST){ + evaluateValues_Const(cub_points,jac_inv,weighted_measure); + } else if(elmtspace == PureBasis::HGRAD){ + evaluateValues_HGrad(cub_points,jac_inv,weighted_measure); + } else if(elmtspace == PureBasis::HCURL){ + evaluateValues_HCurl(cub_points,jac,jac_det,jac_inv,weighted_measure); + } else if(elmtspace == PureBasis::HDIV){ + evaluateValues_HDiv(cub_points,jac,jac_det,weighted_measure); + } else { + TEUCHOS_TEST_FOR_EXCEPT_MSG(true,"panzer::BasisValues2::evaluateValues : Element space not recognized."); + } + + if(use_vertex_coordinates) { + TEUCHOS_TEST_FOR_EXCEPT_MSG(elmtspace == PureBasis::CONST,"panzer::BasisValues2::evaluateValues : Const basis cannot have basis coordinates."); + evaluateBasisCoordinates(vertex_coordinates); + } + +} + +template +void panzer::BasisValues2:: +evaluateValues_Const(const PHX::MDField & cub_points, + const PHX::MDField & jac_inv, + const PHX::MDField & weighted_measure) +{ + + TEUCHOS_ASSERT(getElementSpace() == PureBasis::CONST); + + typedef Intrepid2::FunctionSpaceTools fst; + MDFieldArrayFactory af("",ddims_,true); + + const panzer::PureBasis & basis = *(basis_layout->getBasis()); + + const int num_points = basis_layout->numPoints(); + const int num_basis = basis.cardinality(); + const int num_dim = basis_layout->dimension(); + const int num_cells = basis_layout->numCells(); + + auto cell_basis_scalar = af.buildStaticArray("cell_basis_scalar",1,num_basis,num_points); + auto cell_cub_points = af.buildStaticArray("cell_cub_points",num_points,num_dim); + auto cell_grad_basis = af.buildStaticArray("cell_grad_basis",1,num_basis,num_points,num_dim); + auto cell_jac_inv = af.buildStaticArray("cell_jac_inv",1,num_points,num_dim,num_dim); + + auto cell_basis_ref_scalar = af.buildStaticArray("cell_basis_ref_scalar",num_basis,num_points); + auto cell_grad_basis_ref = af.buildStaticArray("cell_grad_basis_ref",num_basis,num_points,num_dim); + + for(int cell=0;cellgetValues(cell_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE); + + if(compute_derivatives){ + Kokkos::deep_copy(cell_grad_basis_ref.get_view(),0.0); + } + + // ============================================= + // Transform reference values to physical values + + fst::HGRADtransformVALUE(cell_basis_scalar.get_view(),cell_basis_ref_scalar.get_view()); + for(int b=0;b +void panzer::BasisValues2:: +evaluateValues_HGrad(const PHX::MDField & cub_points, + const PHX::MDField & jac_inv, + const PHX::MDField & weighted_measure) +{ + + TEUCHOS_ASSERT(getElementSpace() == PureBasis::HGRAD); + + typedef Intrepid2::FunctionSpaceTools fst; + MDFieldArrayFactory af("",ddims_,true); + + const panzer::PureBasis & basis = *(basis_layout->getBasis()); + + const int num_points = basis_layout->numPoints(); + const int num_basis = basis.cardinality(); + const int num_dim = basis_layout->dimension(); + const int num_cells = cub_points.extent(0); + + auto cell_basis_scalar = af.buildStaticArray("cell_basis_scalar",1,num_basis,num_points); + auto cell_cub_points = af.buildStaticArray("cell_cub_points",num_points,num_dim); + auto cell_grad_basis = af.buildStaticArray("cell_grad_basis",1,num_basis,num_points,num_dim); + auto cell_jac_inv = af.buildStaticArray("cell_jac_inv",1,num_points,num_dim,num_dim); + + auto cell_basis_ref_scalar = af.buildStaticArray("cell_basis_ref_scalar",num_basis,num_points); + auto cell_grad_basis_ref = af.buildStaticArray("cell_grad_basis_ref",num_basis,num_points,num_dim); + + for(int cell=0;cellgetValues(cell_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE); + + if(compute_derivatives){ + intrepid_basis->getValues(cell_grad_basis_ref.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_GRAD); + } + + // ============================================= + // Transform reference values to physical values + + fst::HGRADtransformVALUE(cell_basis_scalar.get_view(),cell_basis_ref_scalar.get_view()); + for(int b=0;b +void panzer::BasisValues2:: +evaluateValues_HCurl(const PHX::MDField & cub_points, + const PHX::MDField & jac, + const PHX::MDField & jac_det, + const PHX::MDField & jac_inv, + const PHX::MDField & weighted_measure) +{ + + TEUCHOS_ASSERT(getElementSpace() == PureBasis::HCURL); + + + typedef Intrepid2::FunctionSpaceTools fst; + MDFieldArrayFactory af("",ddims_,true); + + const panzer::PureBasis & basis = *(basis_layout->getBasis()); + + const int num_points = basis_layout->numPoints(); + const int num_basis = basis.cardinality(); + const int num_dim = basis_layout->dimension(); + const int num_cells = basis_layout->numCells(); + + auto cell_cub_points = af.buildStaticArray("cell_cub_points",num_points,num_dim); + auto cell_jac = af.buildStaticArray("cell_jac",1,num_points,num_dim,num_dim); + auto cell_jac_inv = af.buildStaticArray("cell_jac_inv",1,num_points,num_dim,num_dim); + auto cell_jac_det = af.buildStaticArray("cell_jac_det",1,num_points); + + auto cell_basis_vector = af.buildStaticArray("cell_basis_vector",1,num_basis,num_points,num_dim); + auto cell_curl_basis_scalar = af.buildStaticArray("cell_curl_basis_scalar",1,num_basis,num_points); + auto cell_curl_basis_vector = af.buildStaticArray("cell_curl_basis_vector",1,num_basis,num_points,num_dim); + + auto cell_curl_basis_ref = af.buildArray("cell_curl_basis_ref",num_basis,num_points,num_dim); + auto cell_curl_basis_ref_scalar = af.buildStaticArray("cell_curl_basis_ref_scalar",num_basis,num_points); + auto cell_basis_ref_vector = af.buildArray("cell_basis_ref_vector",num_basis,num_points,num_dim); + + for(int cell=0;cellgetValues(cell_basis_ref_vector.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE); + + if(compute_derivatives){ + if(num_dim==2){ + intrepid_basis->getValues(cell_curl_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_CURL); + } else if(num_dim==3){ + intrepid_basis->getValues(cell_curl_basis_ref.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_CURL); + } + } + + // ============================================= + // Transform reference values to physical values + + fst::HCURLtransformVALUE(cell_basis_vector.get_view(),cell_jac_inv.get_view(),cell_basis_ref_vector.get_view()); + for(int b=0;b +void panzer::BasisValues2:: +evaluateValues_HDiv(const PHX::MDField & cub_points, + const PHX::MDField & jac, + const PHX::MDField & jac_det, + const PHX::MDField & weighted_measure) +{ + + TEUCHOS_ASSERT(getElementSpace() == PureBasis::HDIV); + + typedef Intrepid2::FunctionSpaceTools fst; + MDFieldArrayFactory af("",ddims_,true); + + const panzer::PureBasis & basis = *(basis_layout->getBasis()); + + const int num_points = basis_layout->numPoints(); + const int num_basis = basis.cardinality(); + const int num_dim = basis_layout->dimension(); + const int num_cells = basis_layout->numCells(); + + auto cell_cub_points = af.buildStaticArray("cell_cub_points",num_points,num_dim); + auto cell_jac = af.buildStaticArray("cell_jac",1,num_points,num_dim,num_dim); + auto cell_jac_det = af.buildStaticArray("cell_jac_det",1,num_points); + + auto cell_basis_vector = af.buildStaticArray("cell_basis_vector",1,num_basis,num_points,num_dim); + auto cell_div_basis = af.buildStaticArray("cell_div_basis",1,num_basis,num_points); + + auto cell_basis_ref_vector = af.buildArray("cell_basis_ref_vector",num_basis,num_points,num_dim); + auto cell_div_basis_ref = af.buildStaticArray("cell_div_basis_ref",num_basis,num_points); + + for(int cell=0;cellgetValues(cell_basis_ref_vector.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE); + + if(compute_derivatives){ + intrepid_basis->getValues(cell_div_basis_ref.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_DIV); + } + + // ============================================= + // Transform reference values to physical values + + fst::HDIVtransformVALUE(cell_basis_vector.get_view(),cell_jac.get_view(),cell_jac_det.get_view(),cell_basis_ref_vector.get_view()); + for(int b=0;b +void panzer::BasisValues2:: +evaluateValuesCV(const PHX::MDField & cell_cub_points, + const PHX::MDField & jac, + const PHX::MDField & jac_det, + const PHX::MDField & jac_inv) +{ + MDFieldArrayFactory af("",ddims_,true); + + int num_ip = basis_layout->numPoints(); + int num_card = basis_layout->cardinality(); + int num_dim = basis_layout->dimension(); + + size_type num_cells = jac.extent(0); + + PureBasis::EElementSpace elmtspace = getElementSpace(); + ArrayDynamic dyn_cub_points = af.buildArray("dyn_cub_points", num_ip, num_dim); + + // Integration points are located on physical cells rather than reference cells, + // so we evaluate the basis in a loop over cells. + for (size_type icell = 0; icell < num_cells; ++icell) + { + for (int ip = 0; ip < num_ip; ++ip) + for (int d = 0; d < num_dim; ++d) + dyn_cub_points(ip,d) = cell_cub_points(icell,ip,d); + + if(elmtspace==PureBasis::CONST) { + ArrayDynamic dyn_basis_ref_scalar = af.buildArray("dyn_basis_ref_scalar",num_card,num_ip); + + intrepid_basis->getValues(dyn_basis_ref_scalar.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_VALUE); + + // transform values method just transfers values to array with cell index - no need to call + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_ip; ++ip) + basis_scalar(icell,b,ip) = dyn_basis_ref_scalar(b,ip); + + } + if(elmtspace==PureBasis::HGRAD) { + ArrayDynamic dyn_basis_ref_scalar = af.buildArray("dyn_basis_ref_scalar",num_card,num_ip); + + intrepid_basis->getValues(dyn_basis_ref_scalar.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_VALUE); + + // transform values method just transfers values to array with cell index - no need to call + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_ip; ++ip) + basis_scalar(icell,b,ip) = dyn_basis_ref_scalar(b,ip); + + if(compute_derivatives) { + + int one_cell = 1; + ArrayDynamic dyn_grad_basis_ref = af.buildArray("dyn_grad_basis_ref",num_card,num_ip,num_dim); + ArrayDynamic dyn_grad_basis = af.buildArray("dyn_grad_basis",one_cell,num_card,num_ip,num_dim); + ArrayDynamic dyn_jac_inv = af.buildArray("dyn_jac_inv",one_cell,num_ip,num_dim,num_dim); + + intrepid_basis->getValues(dyn_grad_basis_ref.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_GRAD); + + int cellInd = 0; + for (int ip = 0; ip < num_ip; ++ip) + for (int d1 = 0; d1 < num_dim; ++d1) + for (int d2 = 0; d2 < num_dim; ++d2) + dyn_jac_inv(cellInd,ip,d1,d2) = jac_inv(icell,ip,d1,d2); + + Intrepid2::FunctionSpaceTools::HGRADtransformGRAD(dyn_grad_basis.get_view(), + dyn_jac_inv.get_view(), + dyn_grad_basis_ref.get_view()); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_ip; ++ip) + for (int d = 0; d < num_dim; ++d) + grad_basis(icell,b,ip,d) = dyn_grad_basis(0,b,ip,d); + + } + } + else if(elmtspace==PureBasis::HCURL) { + ArrayDynamic dyn_basis_ref_vector = af.buildArray("dyn_basis_ref_vector",num_card,num_ip,num_dim); + + intrepid_basis->getValues(dyn_basis_ref_vector.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_VALUE); + + int one_cell = 1; + ArrayDynamic dyn_basis_vector = af.buildArray("dyn_basis_vector",one_cell,num_card,num_ip,num_dim); + ArrayDynamic dyn_jac_inv = af.buildArray("dyn_jac_inv",one_cell,num_ip,num_dim,num_dim); + + int cellInd = 0; + for (int ip = 0; ip < num_ip; ++ip) + for (int d1 = 0; d1 < num_dim; ++d1) + for (int d2 = 0; d2 < num_dim; ++d2) + dyn_jac_inv(cellInd,ip,d1,d2) = jac_inv(icell,ip,d1,d2); + + Intrepid2::FunctionSpaceTools::HCURLtransformVALUE(dyn_basis_vector.get_view(), + dyn_jac_inv.get_view(), + dyn_basis_ref_vector.get_view()); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_ip; ++ip) + for (int d = 0; d < num_dim; ++d) + basis_vector(icell,b,ip,d) = dyn_basis_vector(0,b,ip,d); + + if(compute_derivatives && num_dim ==2) { + + int one_cell = 1; + ArrayDynamic dyn_curl_basis_ref_scalar = af.buildArray("dyn_curl_basis_ref_scalar",num_card,num_ip); + ArrayDynamic dyn_curl_basis_scalar = af.buildArray("dyn_curl_basis_scalar",one_cell,num_card,num_ip); + ArrayDynamic dyn_jac_det = af.buildArray("dyn_jac_det",one_cell,num_ip); + + intrepid_basis->getValues(dyn_curl_basis_ref_scalar.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_CURL); + + int cellInd = 0; + for (int ip = 0; ip < num_ip; ++ip) + dyn_jac_det(cellInd,ip) = jac_det(icell,ip); + + Intrepid2::FunctionSpaceTools::HDIVtransformDIV(dyn_curl_basis_scalar.get_view(), + dyn_jac_det.get_view(), + dyn_curl_basis_ref_scalar.get_view()); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_ip; ++ip) + curl_basis_scalar(icell,b,ip) = dyn_curl_basis_scalar(0,b,ip); + + } + if(compute_derivatives && num_dim ==3) { + + int one_cell = 1; + ArrayDynamic dyn_curl_basis_ref = af.buildArray("dyn_curl_basis_ref_vector",num_card,num_ip,num_dim); + ArrayDynamic dyn_curl_basis = af.buildArray("dyn_curl_basis_vector",one_cell,num_card,num_ip,num_dim); + ArrayDynamic dyn_jac_det = af.buildArray("dyn_jac_det",one_cell,num_ip); + ArrayDynamic dyn_jac = af.buildArray("dyn_jac",one_cell,num_ip,num_dim,num_dim); + + intrepid_basis->getValues(dyn_curl_basis_ref.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_CURL); + + int cellInd = 0; + for (int ip = 0; ip < num_ip; ++ip) + { + dyn_jac_det(cellInd,ip) = jac_det(icell,ip); + for (int d1 = 0; d1 < num_dim; ++d1) + for (int d2 = 0; d2 < num_dim; ++d2) + dyn_jac(cellInd,ip,d1,d2) = jac(icell,ip,d1,d2); + } + + Intrepid2::FunctionSpaceTools::HCURLtransformCURL(dyn_curl_basis.get_view(), + dyn_jac.get_view(), + dyn_jac_det.get_view(), + dyn_curl_basis_ref.get_view()); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_ip; ++ip) + for (int d = 0; d < num_dim; ++d) + curl_basis_vector(icell,b,ip,d) = dyn_curl_basis(0,b,ip,d); + + } + + } + else if(elmtspace==PureBasis::HDIV) { + + ArrayDynamic dyn_basis_ref_vector = af.buildArray("dyn_basis_ref_vector",num_card,num_ip,num_dim); + + intrepid_basis->getValues(dyn_basis_ref_vector.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_VALUE); + + int one_cell= 1; + ArrayDynamic dyn_basis_vector = af.buildArray("dyn_basis_vector",one_cell,num_card,num_ip,num_dim); + ArrayDynamic dyn_jac = af.buildArray("dyn_jac",one_cell,num_ip,num_dim,num_dim); + ArrayDynamic dyn_jac_det = af.buildArray("dyn_jac_det",one_cell,num_ip); + + int cellInd = 0; + for (int ip = 0; ip < num_ip; ++ip) + { + dyn_jac_det(cellInd,ip) = jac_det(icell,ip); + for (int d1 = 0; d1 < num_dim; ++d1) + for (int d2 = 0; d2 < num_dim; ++d2) + dyn_jac(cellInd,ip,d1,d2) = jac(icell,ip,d1,d2); + } + + Intrepid2::FunctionSpaceTools::HDIVtransformVALUE(dyn_basis_vector.get_view(), + dyn_jac.get_view(), + dyn_jac_det.get_view(), + dyn_basis_ref_vector.get_view()); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_ip; ++ip) + for (int d = 0; d < num_dim; ++d) + basis_vector(icell,b,ip,d) = dyn_basis_vector(0,b,ip,d); + + if(compute_derivatives) { + + ArrayDynamic dyn_div_basis_ref = af.buildArray("dyn_div_basis_ref_scalar",num_card,num_ip); + ArrayDynamic dyn_div_basis = af.buildArray("dyn_div_basis_scalar",one_cell,num_card,num_ip); + + intrepid_basis->getValues(dyn_div_basis_ref.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_DIV); + + Intrepid2::FunctionSpaceTools::HDIVtransformDIV(dyn_div_basis.get_view(), + dyn_jac_det.get_view(), + dyn_div_basis_ref.get_view()); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_ip; ++ip) + div_basis(icell,b,ip) = dyn_div_basis(0,b,ip); + + } + + } + else { TEUCHOS_ASSERT(false); } + + } // cell loop + +} + +template +void panzer::BasisValues2:: +evaluateReferenceValues(const PHX::MDField & cub_points,bool compute_derivatives,bool use_vertex_coordinates) +{ + MDFieldArrayFactory af("",ddims_,true); + + int num_quad = basis_layout->numPoints(); + int num_dim = basis_layout->dimension(); + int num_card = basis_layout->cardinality(); + + ArrayDynamic dyn_cub_points = af.buildArray("dyn_cub_points", num_quad,num_dim); + + for (int ip = 0; ip < num_quad; ++ip) + for (int d = 0; d < num_dim; ++d) + dyn_cub_points(ip,d) = cub_points(ip,d); + + PureBasis::EElementSpace elmtspace = getElementSpace(); + if(elmtspace==PureBasis::HGRAD || elmtspace==PureBasis::CONST) { + ArrayDynamic dyn_basis_ref_scalar = af.buildArray("dyn_basis_ref_scalar",num_card,num_quad); + + intrepid_basis->getValues(dyn_basis_ref_scalar.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_VALUE); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_quad; ++ip) + basis_ref_scalar(b,ip) = dyn_basis_ref_scalar(b,ip); + } + else if(elmtspace==PureBasis::HDIV || elmtspace==PureBasis::HCURL) { + ArrayDynamic dyn_basis_ref_vector = af.buildArray("dyn_basis_ref_vector",num_card,num_quad,num_dim); + + intrepid_basis->getValues(dyn_basis_ref_vector.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_VALUE); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_quad; ++ip) + for (int d = 0; d < num_dim; ++d) + basis_ref_vector(b,ip,d) = dyn_basis_ref_vector(b,ip,d); + } + else { TEUCHOS_ASSERT(false); } + + if(elmtspace==PureBasis::HGRAD && compute_derivatives) { + ArrayDynamic dyn_grad_basis_ref = af.buildArray("dyn_basis_ref_vector",num_card,num_quad,num_dim); + + intrepid_basis->getValues(dyn_grad_basis_ref.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_GRAD); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_quad; ++ip) + for (int d = 0; d < num_dim; ++d) + grad_basis_ref(b,ip,d) = dyn_grad_basis_ref(b,ip,d); + } + else if(elmtspace==PureBasis::HCURL && compute_derivatives && num_dim==2) { + ArrayDynamic dyn_curl_basis_ref = af.buildArray("dyn_curl_basis_ref_scalar",num_card,num_quad); + + intrepid_basis->getValues(dyn_curl_basis_ref.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_CURL); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_quad; ++ip) + curl_basis_ref_scalar(b,ip) = dyn_curl_basis_ref(b,ip); + } + else if(elmtspace==PureBasis::HCURL && compute_derivatives && num_dim==3) { + ArrayDynamic dyn_curl_basis_ref = af.buildArray("dyn_curl_basis_ref_vector",num_card,num_quad,num_dim); + + intrepid_basis->getValues(dyn_curl_basis_ref.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_CURL); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_quad; ++ip) + for (int d = 0; d < num_dim; ++d) + curl_basis_ref_vector(b,ip,d) = dyn_curl_basis_ref(b,ip,d); + } + else if(elmtspace==PureBasis::HDIV && compute_derivatives) { + ArrayDynamic dyn_div_basis_ref = af.buildArray("dyn_div_basis_ref_scalar",num_card,num_quad); + + intrepid_basis->getValues(dyn_div_basis_ref.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_DIV); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_quad; ++ip) + div_basis_ref(b,ip) = dyn_div_basis_ref(b,ip); + } + + + if(use_vertex_coordinates) { + // Intrepid removes fad types from the coordinate scalar type. We + // pull the actual field scalar type from the basis object to be + // consistent. + if (elmtspace != PureBasis::CONST) { + using coordsScalarType = typename Intrepid2::Basis::scalarType; + auto dyn_basis_coordinates_ref = af.buildArray("basis_coordinates_ref", + basis_coordinates_ref.extent(0), + basis_coordinates_ref.extent(1)); + intrepid_basis->getDofCoords(dyn_basis_coordinates_ref.get_view()); + + // fill in basis coordinates + for (int i = 0; i < basis_coordinates_ref.extent_int(0); ++i) + for (int j = 0; j < basis_coordinates_ref.extent_int(1); ++j) + basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j); + } + } + + references_evaluated = true; +} + +// method for applying orientations +template +void BasisValues2:: +applyOrientations(const std::vector & orientations) +{ + if (!intrepid_basis->requireOrientation()) + return; + + typedef Intrepid2::OrientationTools ots; + const PureBasis::EElementSpace elmtspace = getElementSpace(); + + // orientation (right now std vector) is created using push_back method. + // thus, its size is the actual number of elements to be applied. + // on the other hand, basis_layout num cell indicates workset size. + // to get right size of cells, use minimum of them. + const int num_cell_basis_layout = basis_layout->numCells(), num_cell_orientation = orientations.size(); + const int num_cell = num_cell_basis_layout < num_cell_orientation ? num_cell_basis_layout : num_cell_orientation; + const int num_dim = basis_layout->dimension(); + const Kokkos::pair range_cell(0, num_cell); + + // Kokkos::DynRankView + // drv_orts((Intrepid2::Orientation*)orientations.data(), num_cell); + Kokkos::DynRankView drv_orts("drv_orts", num_cell); + auto host_drv_orts = Kokkos::create_mirror_view(drv_orts); + for (size_t i=0; i < drv_orts.size(); ++i) + host_drv_orts(i) = orientations[i]; + Kokkos::deep_copy(drv_orts,host_drv_orts); + PHX::Device::fence(); + + /// + /// HGRAD elements + /// + if (elmtspace==PureBasis::HGRAD) { + { + { + auto drv_basis_scalar = Kokkos::subview(basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); + auto drv_basis_scalar_tmp = Kokkos::createDynRankView(basis_scalar.get_view(), + "drv_basis_scalar_tmp", + drv_basis_scalar.extent(0), // C + drv_basis_scalar.extent(1), // F + drv_basis_scalar.extent(2)); // P + Kokkos::deep_copy(drv_basis_scalar_tmp, drv_basis_scalar); + ots::modifyBasisByOrientation(drv_basis_scalar, + drv_basis_scalar_tmp, + drv_orts, + intrepid_basis); + } + if(build_weighted) { + auto drv_basis_scalar = Kokkos::subview(weighted_basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); + auto drv_basis_scalar_tmp = Kokkos::createDynRankView(weighted_basis_scalar.get_view(), + "drv_basis_scalar_tmp", + drv_basis_scalar.extent(0), // C + drv_basis_scalar.extent(1), // F + drv_basis_scalar.extent(2)); // P + Kokkos::deep_copy(drv_basis_scalar_tmp, drv_basis_scalar); + ots::modifyBasisByOrientation(drv_basis_scalar, + drv_basis_scalar_tmp, + drv_orts, + intrepid_basis); + } + + } + + if (compute_derivatives) { + { + auto drv_grad_basis = Kokkos::subview(grad_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + auto drv_grad_basis_tmp = Kokkos::createDynRankView(grad_basis.get_view(), + "drv_grad_basis_tmp", + drv_grad_basis.extent(0), // C + drv_grad_basis.extent(1), // F + drv_grad_basis.extent(2), // P + drv_grad_basis.extent(3)); // D + Kokkos::deep_copy(drv_grad_basis_tmp, drv_grad_basis); + ots::modifyBasisByOrientation(drv_grad_basis, + drv_grad_basis_tmp, + drv_orts, + intrepid_basis); + } + if(build_weighted) { + auto drv_grad_basis = Kokkos::subview(weighted_grad_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + auto drv_grad_basis_tmp = Kokkos::createDynRankView(weighted_grad_basis.get_view(), + "drv_grad_basis_tmp", + drv_grad_basis.extent(0), // C + drv_grad_basis.extent(1), // F + drv_grad_basis.extent(2), // P + drv_grad_basis.extent(3)); // D + Kokkos::deep_copy(drv_grad_basis_tmp, drv_grad_basis); + ots::modifyBasisByOrientation(drv_grad_basis, + drv_grad_basis_tmp, + drv_orts, + intrepid_basis); + } + } + } + + /// + /// hcurl 2d elements + /// + else if (elmtspace==PureBasis::HCURL && num_dim==2) { + { + { + auto drv_basis_vector = Kokkos::subview(basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + auto drv_basis_vector_tmp = Kokkos::createDynRankView(basis_vector.get_view(), + "drv_basis_vector_tmp", + drv_basis_vector.extent(0), // C + drv_basis_vector.extent(1), // F + drv_basis_vector.extent(2), // P + drv_basis_vector.extent(3)); // D + Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); + ots::modifyBasisByOrientation(drv_basis_vector, + drv_basis_vector_tmp, + drv_orts, + intrepid_basis); + } + if(build_weighted) { + auto drv_basis_vector = Kokkos::subview(weighted_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + auto drv_basis_vector_tmp = Kokkos::createDynRankView(weighted_basis_vector.get_view(), + "drv_basis_vector_tmp", + drv_basis_vector.extent(0), // C + drv_basis_vector.extent(1), // F + drv_basis_vector.extent(2), // P + drv_basis_vector.extent(3)); // D + Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); + ots::modifyBasisByOrientation(drv_basis_vector, + drv_basis_vector_tmp, + drv_orts, + intrepid_basis); + } + } + + if (compute_derivatives) { + { + auto drv_curl_basis_scalar = Kokkos::subview(curl_basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); + auto drv_curl_basis_scalar_tmp = Kokkos::createDynRankView(curl_basis_scalar.get_view(), + "drv_curl_basis_scalar_tmp", + drv_curl_basis_scalar.extent(0), // C + drv_curl_basis_scalar.extent(1), // F + drv_curl_basis_scalar.extent(2)); // P + Kokkos::deep_copy(drv_curl_basis_scalar_tmp, drv_curl_basis_scalar); + ots::modifyBasisByOrientation(drv_curl_basis_scalar, + drv_curl_basis_scalar_tmp, + drv_orts, + intrepid_basis); + } + + if(build_weighted) { + auto drv_curl_basis_scalar = Kokkos::subview(weighted_curl_basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); + auto drv_curl_basis_scalar_tmp = Kokkos::createDynRankView(weighted_curl_basis_scalar.get_view(), + "drv_curl_basis_scalar_tmp", + drv_curl_basis_scalar.extent(0), // C + drv_curl_basis_scalar.extent(1), // F + drv_curl_basis_scalar.extent(2)); // P + Kokkos::deep_copy(drv_curl_basis_scalar_tmp, drv_curl_basis_scalar); + ots::modifyBasisByOrientation(drv_curl_basis_scalar, + drv_curl_basis_scalar_tmp, + drv_orts, + intrepid_basis); + } + } + } + + /// + /// hcurl 3d elements + /// + else if (elmtspace==PureBasis::HCURL && num_dim==3) { + { + { + auto drv_basis_vector = Kokkos::subview(basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + auto drv_basis_vector_tmp = Kokkos::createDynRankView(basis_vector.get_view(), + "drv_basis_vector_tmp", + drv_basis_vector.extent(0), // C + drv_basis_vector.extent(1), // F + drv_basis_vector.extent(2), // P + drv_basis_vector.extent(3)); // D + Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); + ots::modifyBasisByOrientation(drv_basis_vector, + drv_basis_vector_tmp, + drv_orts, + intrepid_basis); + } + if(build_weighted) { + auto drv_basis_vector = Kokkos::subview(weighted_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + auto drv_basis_vector_tmp = Kokkos::createDynRankView(weighted_basis_vector.get_view(), + "drv_basis_vector_tmp", + drv_basis_vector.extent(0), // C + drv_basis_vector.extent(1), // F + drv_basis_vector.extent(2), // P + drv_basis_vector.extent(3)); // D + Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); + ots::modifyBasisByOrientation(drv_basis_vector, + drv_basis_vector_tmp, + drv_orts, + intrepid_basis); + } + } + + if (compute_derivatives) { + { + auto drv_curl_basis_vector = Kokkos::subview(curl_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + auto drv_curl_basis_vector_tmp = Kokkos::createDynRankView(curl_basis_vector.get_view(), + "drv_curl_basis_vector_tmp", + drv_curl_basis_vector.extent(0), // C + drv_curl_basis_vector.extent(1), // F + drv_curl_basis_vector.extent(2), // P + drv_curl_basis_vector.extent(3)); // D + Kokkos::deep_copy(drv_curl_basis_vector_tmp, drv_curl_basis_vector); + ots::modifyBasisByOrientation(drv_curl_basis_vector, + drv_curl_basis_vector_tmp, + drv_orts, + intrepid_basis); + } + if(build_weighted) { + auto drv_curl_basis_vector = Kokkos::subview(weighted_curl_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + auto drv_curl_basis_vector_tmp = Kokkos::createDynRankView(weighted_curl_basis_vector.get_view(), + "drv_curl_basis_vector_tmp", + drv_curl_basis_vector.extent(0), // C + drv_curl_basis_vector.extent(1), // F + drv_curl_basis_vector.extent(2), // P + drv_curl_basis_vector.extent(3)); // D + Kokkos::deep_copy(drv_curl_basis_vector_tmp, drv_curl_basis_vector); + ots::modifyBasisByOrientation(drv_curl_basis_vector, + drv_curl_basis_vector_tmp, + drv_orts, + intrepid_basis); + } + } + } + /// + /// hdiv elements (2d and 3d) + /// + else if (elmtspace==PureBasis::HDIV) { + { + { + auto drv_basis_vector = Kokkos::subview(basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + auto drv_basis_vector_tmp = Kokkos::createDynRankView(basis_vector.get_view(), + "drv_basis_vector_tmp", + drv_basis_vector.extent(0), // C + drv_basis_vector.extent(1), // F + drv_basis_vector.extent(2), // P + drv_basis_vector.extent(3)); // D + Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); + ots::modifyBasisByOrientation(drv_basis_vector, + drv_basis_vector_tmp, + drv_orts, + intrepid_basis); + } + if(build_weighted) { + auto drv_basis_vector = Kokkos::subview(weighted_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + auto drv_basis_vector_tmp = Kokkos::createDynRankView(weighted_basis_vector.get_view(), + "drv_basis_vector_tmp", + drv_basis_vector.extent(0), // C + drv_basis_vector.extent(1), // F + drv_basis_vector.extent(2), // P + drv_basis_vector.extent(3)); // D + Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector); + ots::modifyBasisByOrientation(drv_basis_vector, + drv_basis_vector_tmp, + drv_orts, + intrepid_basis); + } + } + if (compute_derivatives) { + { + auto drv_div_basis = Kokkos::subview(div_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); + auto drv_div_basis_tmp = Kokkos::createDynRankView(div_basis.get_view(), + "drv_div_basis_tmp", + drv_div_basis.extent(0), // C + drv_div_basis.extent(1), // F + drv_div_basis.extent(2)); // P + Kokkos::deep_copy(drv_div_basis_tmp, drv_div_basis); + ots::modifyBasisByOrientation(drv_div_basis, + drv_div_basis_tmp, + drv_orts, + intrepid_basis); + } + if(build_weighted) { + auto drv_div_basis = Kokkos::subview(weighted_div_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL()); + auto drv_div_basis_tmp = Kokkos::createDynRankView(weighted_div_basis.get_view(), + "drv_div_basis_tmp", + drv_div_basis.extent(0), // C + drv_div_basis.extent(1), // F + drv_div_basis.extent(2)); // P + Kokkos::deep_copy(drv_div_basis_tmp, drv_div_basis); + ots::modifyBasisByOrientation(drv_div_basis, + drv_div_basis_tmp, + drv_orts, + intrepid_basis); + } + } + } +} + +// method for applying orientations +template +void BasisValues2:: +applyOrientations(const PHX::MDField & orientations) +{ + + TEUCHOS_TEST_FOR_EXCEPT_MSG(true,"panzer::BasisValues2::applyOrientations : this should not be called."); + + int num_cell = orientations.extent(0); + int num_basis = orientations.extent(1); + int num_dim = basis_layout->dimension(); + int num_ip = basis_layout->numPoints(); + PureBasis::EElementSpace elmtspace = getElementSpace(); + + if(elmtspace==PureBasis::HCURL && num_dim==2) { + + // setup the orientations for the trial space + // Intrepid2::FunctionSpaceTools::applyFieldSigns(basis_vector,orientations); + + for (int c=0; c(curl_basis_scalar,orientations); + for (int c=0; c::applyFieldSigns(weighted_basis_vector.get_view(),orientations.get_view()); + + if(compute_derivatives) + Intrepid2::FunctionSpaceTools::applyFieldSigns(weighted_curl_basis_scalar.get_view(),orientations.get_view()); + } + } + else if(elmtspace==PureBasis::HCURL && num_dim==3) { + + // setup the orientations for the trial space + // Intrepid2::FunctionSpaceTools::applyFieldSigns(basis_vector,orientations); + + for (int c=0; c(curl_basis_vector,orientations); + for (int c=0; c::applyFieldSigns(weighted_basis_vector.get_view(),orientations.get_view()); + + if(compute_derivatives) + Intrepid2::FunctionSpaceTools::applyFieldSigns(weighted_curl_basis_vector.get_view(),orientations.get_view()); + } + } + else if(elmtspace==PureBasis::HDIV) { + // setup the orientations for the trial space + // Intrepid2::FunctionSpaceTools::applyFieldSigns(basis_vector,orientations); + + for (int c=0; c(div_basis,orientations); + + for (int c=0; c::applyFieldSigns(weighted_basis_vector.get_view(),orientations.get_view()); + + if(compute_derivatives) + Intrepid2::FunctionSpaceTools::applyFieldSigns(weighted_div_basis.get_view(),orientations.get_view()); + } + } +} + +template +PureBasis::EElementSpace BasisValues2::getElementSpace() const +{ return basis_layout->getBasis()->getElementSpace(); } + +template +void panzer::BasisValues2:: +setupArrays(const Teuchos::RCP& layout, + bool computeDerivatives) +{ + MDFieldArrayFactory af(prefix,alloc_arrays); + + compute_derivatives = computeDerivatives; + basis_layout = layout; + Teuchos::RCP basisDesc = layout->getBasis(); + + // for convience pull out basis and quadrature information + int num_quad = layout->numPoints(); + int dim = basisDesc->dimension(); + int card = basisDesc->cardinality(); + int numcells = basisDesc->numCells(); + panzer::PureBasis::EElementSpace elmtspace = basisDesc->getElementSpace(); + Teuchos::RCP cellTopo = basisDesc->getCellTopology(); + + intrepid_basis = basisDesc->getIntrepid2Basis(); + + // allocate field containers + // field sizes defined by http://trilinos.sandia.gov/packages/docs/dev/packages/intrepid/doc/html/basis_page.html#basis_md_array_sec + + // compute basis fields + if(elmtspace==panzer::PureBasis::HGRAD) { + // HGRAD is a nodal field + + // build values + /////////////////////////////////////////////// + basis_ref_scalar = af.buildStaticArray("basis_ref",card,num_quad); // F, P + basis_scalar = af.buildStaticArray("basis",numcells,card,num_quad); + + if(build_weighted) + weighted_basis_scalar = af.buildStaticArray("weighted_basis",numcells,card,num_quad); + + // build gradients + /////////////////////////////////////////////// + + if(compute_derivatives) { + grad_basis_ref = af.buildStaticArray("grad_basis_ref",card,num_quad,dim); // F, P, D + grad_basis = af.buildStaticArray("grad_basis",numcells,card,num_quad,dim); + + if(build_weighted) + weighted_grad_basis = af.buildStaticArray("weighted_grad_basis",numcells,card,num_quad,dim); + } + + // build curl + /////////////////////////////////////////////// + + // nothing - HGRAD does not support CURL operation + } + else if(elmtspace==panzer::PureBasis::HCURL) { + // HCURL is a vector field + + // build values + /////////////////////////////////////////////// + + basis_ref_vector = af.buildStaticArray("basis_ref",card,num_quad,dim); // F, P, D + basis_vector = af.buildStaticArray("basis",numcells,card,num_quad,dim); + + if(build_weighted) + weighted_basis_vector = af.buildStaticArray("weighted_basis",numcells,card,num_quad,dim); + + // build gradients + /////////////////////////////////////////////// + + // nothing - HCURL does not support GRAD operation + + // build curl + /////////////////////////////////////////////// + + if(compute_derivatives) { + if(dim==2) { + // curl of HCURL basis is not dimension dependent + curl_basis_ref_scalar = af.buildStaticArray("curl_basis_ref",card,num_quad); // F, P + curl_basis_scalar = af.buildStaticArray("curl_basis",numcells,card,num_quad); + + if(build_weighted) + weighted_curl_basis_scalar = af.buildStaticArray("weighted_curl_basis",numcells,card,num_quad); + } + else if(dim==3){ + curl_basis_ref_vector = af.buildStaticArray("curl_basis_ref",card,num_quad,dim); // F, P, D + curl_basis_vector = af.buildStaticArray("curl_basis",numcells,card,num_quad,dim); + + if(build_weighted) + weighted_curl_basis_vector = af.buildStaticArray("weighted_curl_basis",numcells,card,num_quad,dim); + } + else { TEUCHOS_ASSERT(false); } // what do I do with 1D? + } + } + else if(elmtspace==panzer::PureBasis::HDIV) { + // HDIV is a vector field + + // build values + /////////////////////////////////////////////// + + basis_ref_vector = af.buildStaticArray("basis_ref",card,num_quad,dim); // F, P, D + basis_vector = af.buildStaticArray("basis",numcells,card,num_quad,dim); + + if(build_weighted) + weighted_basis_vector = af.buildStaticArray("weighted_basis",numcells,card,num_quad,dim); + + // build gradients + /////////////////////////////////////////////// + + // nothing - HCURL does not support GRAD operation + + // build curl + /////////////////////////////////////////////// + + // nothing - HDIV does not support CURL operation + + // build div + /////////////////////////////////////////////// + + if(compute_derivatives) { + div_basis_ref = af.buildStaticArray("div_basis_ref",card,num_quad); // F, P + div_basis = af.buildStaticArray("div_basis",numcells,card,num_quad); + + if(build_weighted) + weighted_div_basis = af.buildStaticArray("weighted_div_basis",numcells,card,num_quad); + } + } + else if(elmtspace==panzer::PureBasis::CONST) { + // CONST is a nodal field + + // build values + /////////////////////////////////////////////// + basis_ref_scalar = af.buildStaticArray("basis_ref",card,num_quad); // F, P + basis_scalar = af.buildStaticArray("basis",numcells,card,num_quad); + + if(build_weighted) + weighted_basis_scalar = af.buildStaticArray("weighted_basis",numcells,card,num_quad); + + // build gradients + /////////////////////////////////////////////// + + // nothing - CONST does not support GRAD operation + + // build curl + /////////////////////////////////////////////// + + // nothing - CONST does not support CURL operation + + // build div + /////////////////////////////////////////////// + + // nothing - CONST does not support DIV operation + } + else { TEUCHOS_ASSERT(false); } + + basis_coordinates_ref = af.buildStaticArray("basis_coordinates_ref",card,dim); + basis_coordinates = af.buildStaticArray("basis_coordinates",numcells,card,dim); +} + +} // namespace panzer diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_BasisValues_Evaluator_decl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_BasisValues_Evaluator_decl.hpp index bd731f5c4d71..5d7b11cd83f9 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_BasisValues_Evaluator_decl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_BasisValues_Evaluator_decl.hpp @@ -81,9 +81,9 @@ class BasisValues_Evaluator Teuchos::RCP basis; // is anything other than ScalarT really needed here? - Teuchos::RCP > basisValues; - PointValues2 pointValues; - PointValues2 constPointValues; + Teuchos::RCP > basisValues; + PointValues2 pointValues; + PointValues2 constPointValues; Teuchos::RCP > orientations; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_BasisValues_Evaluator_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_BasisValues_Evaluator_impl.hpp index c185d84a2003..8e4e83c24810 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_BasisValues_Evaluator_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_BasisValues_Evaluator_impl.hpp @@ -100,7 +100,7 @@ void BasisValues_Evaluator::initialize(const Teuchos::RCPdimension(); // setup all fields to be evaluated and constructed - pointValues = PointValues2(pointRule->getName()+"_",false); + pointValues = PointValues2(pointRule->getName()+"_",false); pointValues.setupArrays(pointRule); // the field manager will allocate all of these field @@ -114,7 +114,7 @@ void BasisValues_Evaluator::initialize(const Teuchos::RCP layout = Teuchos::rcp(new panzer::BasisIRLayout(basis,*pointRule)); - basisValues = Teuchos::rcp(new BasisValues2(basis->name()+"_"+pointRule->getName()+"_",false)); + basisValues = Teuchos::rcp(new BasisValues2(basis->name()+"_"+pointRule->getName()+"_",false)); basisValues->setupArrays(layout,derivativesRequired_); // the field manager will allocate all of these field diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_DOF_PointValues.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_DOF_PointValues.hpp index f5d2b7878dcf..8210aa56e4d7 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_DOF_PointValues.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_DOF_PointValues.hpp @@ -76,14 +76,14 @@ class DOF_PointValues : public panzer::EvaluatorWithBaseImpl, bool is_vector_basis; Teuchos::RCP basis; - Teuchos::RCP > basisValues; - PHX::MDField + Teuchos::RCP > basisValues; + PHX::MDField constBasisRefScalar_; - PHX::MDField + PHX::MDField constBasisScalar_; - PHX::MDField + PHX::MDField constBasisRefVector_; - PHX::MDField + PHX::MDField constBasisVector_; }; @@ -117,14 +117,14 @@ class DOF_PointValues : Kokkos::View offsets_array; Teuchos::RCP basis; - Teuchos::RCP > basisValues; - PHX::MDField + Teuchos::RCP > basisValues; + PHX::MDField constBasisRefScalar_; - PHX::MDField + PHX::MDField constBasisScalar_; - PHX::MDField + PHX::MDField constBasisRefVector_; - PHX::MDField + PHX::MDField constBasisVector_; }; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_DOF_PointValues_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_DOF_PointValues_impl.hpp index 0c27bb29853d..ffb9418874f5 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_DOF_PointValues_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_DOF_PointValues_impl.hpp @@ -84,7 +84,7 @@ DOF_PointValues(const Teuchos::ParameterList & p) // setup all basis fields that are required Teuchos::RCP layout = Teuchos::rcp(new BasisIRLayout(basis,*pointRule)); - basisValues = Teuchos::rcp(new BasisValues2(basis->name()+"_"+pointRule->getName()+"_")); + basisValues = Teuchos::rcp(new BasisValues2(basis->name()+"_"+pointRule->getName()+"_")); basisValues->setupArrays(layout,false); // the field manager will allocate all of these field @@ -140,16 +140,16 @@ evaluateFields(typename TRAITS::EvalData workset) if(is_vector_basis) { int spaceDim = basisValues->basis_vector.extent(3); if(spaceDim==3) { - dof_functors::EvaluateDOFWithSens_Vector::Array_CellBasisIPDim,3> functor(dof_basis,dof_ip_vector,basisValues->basis_vector); + dof_functors::EvaluateDOFWithSens_Vector::Array_CellBasisIPDim,3> functor(dof_basis,dof_ip_vector,basisValues->basis_vector); Kokkos::parallel_for(workset.num_cells,functor); } else { - dof_functors::EvaluateDOFWithSens_Vector::Array_CellBasisIPDim,2> functor(dof_basis,dof_ip_vector,basisValues->basis_vector); + dof_functors::EvaluateDOFWithSens_Vector::Array_CellBasisIPDim,2> functor(dof_basis,dof_ip_vector,basisValues->basis_vector); Kokkos::parallel_for(workset.num_cells,functor); } } else { - dof_functors::EvaluateDOFWithSens_Scalar::Array_CellBasisIP> functor(dof_basis,dof_ip_scalar,basisValues->basis_scalar); + dof_functors::EvaluateDOFWithSens_Scalar::Array_CellBasisIP> functor(dof_basis,dof_ip_scalar,basisValues->basis_scalar); Kokkos::parallel_for(workset.num_cells,functor); } } @@ -195,7 +195,7 @@ DOF_PointValues(const Teuchos::ParameterList & p) // setup all basis fields that are required Teuchos::RCP layout = Teuchos::rcp(new BasisIRLayout(basis,*pointRule)); - basisValues = Teuchos::rcp(new BasisValues2(basis->name()+"_"+pointRule->getName()+"_")); + basisValues = Teuchos::rcp(new BasisValues2(basis->name()+"_"+pointRule->getName()+"_")); basisValues->setupArrays(layout,false); // the field manager will allocate all of these field @@ -252,33 +252,33 @@ evaluateFields(typename TRAITS::EvalData workset) if(accelerate_jacobian) { int spaceDim = basisValues->basis_vector.extent(3); if(spaceDim==3) { - dof_functors::EvaluateDOFFastSens_Vector::Array_CellBasisIPDim,3> functor(dof_basis,dof_ip_vector,offsets_array,basisValues->basis_vector); + dof_functors::EvaluateDOFFastSens_Vector::Array_CellBasisIPDim,3> functor(dof_basis,dof_ip_vector,offsets_array,basisValues->basis_vector); Kokkos::parallel_for(workset.num_cells,functor); } else { - dof_functors::EvaluateDOFFastSens_Vector::Array_CellBasisIPDim,2> functor(dof_basis,dof_ip_vector,offsets_array,basisValues->basis_vector); + dof_functors::EvaluateDOFFastSens_Vector::Array_CellBasisIPDim,2> functor(dof_basis,dof_ip_vector,offsets_array,basisValues->basis_vector); Kokkos::parallel_for(workset.num_cells,functor); } } else { int spaceDim = basisValues->basis_vector.extent(3); if(spaceDim==3) { - dof_functors::EvaluateDOFWithSens_Vector::Array_CellBasisIPDim,3> functor(dof_basis,dof_ip_vector,basisValues->basis_vector); + dof_functors::EvaluateDOFWithSens_Vector::Array_CellBasisIPDim,3> functor(dof_basis,dof_ip_vector,basisValues->basis_vector); Kokkos::parallel_for(workset.num_cells,functor); } else { - dof_functors::EvaluateDOFWithSens_Vector::Array_CellBasisIPDim,2> functor(dof_basis,dof_ip_vector,basisValues->basis_vector); + dof_functors::EvaluateDOFWithSens_Vector::Array_CellBasisIPDim,2> functor(dof_basis,dof_ip_vector,basisValues->basis_vector); Kokkos::parallel_for(workset.num_cells,functor); } } } else { if(accelerate_jacobian) { - dof_functors::EvaluateDOFFastSens_Scalar::Array_CellBasisIP> functor(dof_basis,dof_ip_scalar,offsets_array,basisValues->basis_scalar); + dof_functors::EvaluateDOFFastSens_Scalar::Array_CellBasisIP> functor(dof_basis,dof_ip_scalar,offsets_array,basisValues->basis_scalar); Kokkos::parallel_for(workset.num_cells,functor); } else { - dof_functors::EvaluateDOFWithSens_Scalar::Array_CellBasisIP> functor(dof_basis,dof_ip_scalar,basisValues->basis_scalar); + dof_functors::EvaluateDOFWithSens_Scalar::Array_CellBasisIP> functor(dof_basis,dof_ip_scalar,basisValues->basis_scalar); Kokkos::parallel_for(workset.num_cells,functor); } } diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_EdgeBasis.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_EdgeBasis.hpp index 61569af8c2d3..2760a9699fa2 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_EdgeBasis.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_EdgeBasis.hpp @@ -95,8 +95,8 @@ class DirichletResidual_EdgeBasis Teuchos::RCP basis; Teuchos::RCP pointRule; - PointValues2 pointValues; - PHX::MDField + PointValues2 pointValues; + PHX::MDField constJac_; Teuchos::RCP > orientations; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_EdgeBasis_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_EdgeBasis_impl.hpp index 231fedae792b..97fe7d98b6b5 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_EdgeBasis_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_EdgeBasis_impl.hpp @@ -93,7 +93,7 @@ DirichletResidual_EdgeBasis( // setup all basis fields that are required // setup all fields to be evaluated and constructed - pointValues = PointValues2(pointRule->getName()+"_",false); + pointValues = PointValues2(pointRule->getName()+"_",false); pointValues.setupArrays(pointRule); // the field manager will allocate all of these field diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_FaceBasis.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_FaceBasis.hpp index fb8128ef55ef..0d9033a19dc5 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_FaceBasis.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_FaceBasis.hpp @@ -97,8 +97,8 @@ class DirichletResidual_FaceBasis Kokkos::DynRankView faceNormal; // face normals Kokkos::DynRankView refFaceNormal; // reference face normals - PointValues2 pointValues; - PHX::MDField + PointValues2 pointValues; + PHX::MDField constJac_; Teuchos::RCP > orientations; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_FaceBasis_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_FaceBasis_impl.hpp index d333362ec297..65b301d38b5a 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_FaceBasis_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_Dirichlet_Residual_FaceBasis_impl.hpp @@ -86,7 +86,7 @@ DirichletResidual_FaceBasis( value = PHX::MDField(value_name, vector_layout_vector); // setup all fields to be evaluated and constructed - pointValues = PointValues2 (pointRule->getName()+"_",false); + pointValues = PointValues2 (pointRule->getName()+"_",false); pointValues.setupArrays(pointRule); // the field manager will allocate all of these field diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherNormals_decl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherNormals_decl.hpp index ca8946ee05ed..580857eae7f9 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherNormals_decl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherNormals_decl.hpp @@ -93,8 +93,8 @@ class GatherNormals Kokkos::DynRankView faceNormal; // face normals Kokkos::DynRankView refFaceNormal; // reference face normals - PointValues2 pointValues; - PHX::MDField + PointValues2 pointValues; + PHX::MDField constJac_; Teuchos::RCP > orientations; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherNormals_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherNormals_impl.hpp index 89b80366299c..5e8c9e45a5b4 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherNormals_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherNormals_impl.hpp @@ -78,7 +78,7 @@ GatherNormals( // setup the orientation field std::string orientationFieldName = basis->name() + " Orientation"; // setup all fields to be evaluated and constructed - pointValues = panzer::PointValues2 (pointRule->getName()+"_",false); + pointValues = panzer::PointValues2 (pointRule->getName()+"_",false); pointValues.setupArrays(pointRule); // the field manager will allocate all of these field diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangents_decl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangents_decl.hpp index 0e0d88eba333..e3eea502e733 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangents_decl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangents_decl.hpp @@ -93,8 +93,8 @@ class GatherTangents Kokkos::DynRankView edgeTan; // edge tangents Kokkos::DynRankView refEdgeTan; // reference edge tangents - PointValues2 pointValues; - PHX::MDField + PointValues2 pointValues; + PHX::MDField constJac_; Teuchos::RCP > orientations; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangents_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangents_impl.hpp index 5fddf5991122..21346dae0a29 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangents_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangents_impl.hpp @@ -78,7 +78,7 @@ GatherTangents( // setup the orientation field std::string orientationFieldName = basis->name() + " Orientation"; // setup all fields to be evaluated and constructed - pointValues = panzer::PointValues2 (pointRule->getName()+"_",false); + pointValues = panzer::PointValues2 (pointRule->getName()+"_",false); pointValues.setupArrays(pointRule); // the field manager will allocate all of these field diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_PointValues_Evaluator_decl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_PointValues_Evaluator_decl.hpp index b1331f4bf2ed..7151807f5801 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_PointValues_Evaluator_decl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_PointValues_Evaluator_decl.hpp @@ -77,7 +77,7 @@ class PointValues_Evaluator using ScalarT = typename EvalT::ScalarT; // is anything other than ScalarT really needed here? - PointValues2 pointValues; + PointValues2 pointValues; PHX::MDField refPointArray; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_PointValues_Evaluator_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_PointValues_Evaluator_impl.hpp index 4c9fb5c65e52..00ff4fc40099 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_PointValues_Evaluator_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_PointValues_Evaluator_impl.hpp @@ -128,7 +128,7 @@ void PointValues_Evaluator::initialize(const Teuchos::RCP(pointRule->getName()+"_",false); + pointValues = PointValues2(pointRule->getName()+"_",false); pointValues.setupArrays(pointRule); // the field manager will allocate all of these field diff --git a/packages/panzer/disc-fe/test/core_tests/basis_values2.cpp b/packages/panzer/disc-fe/test/core_tests/basis_values2.cpp index c52357611b6a..910cc7dfb41e 100644 --- a/packages/panzer/disc-fe/test/core_tests/basis_values2.cpp +++ b/packages/panzer/disc-fe/test/core_tests/basis_values2.cpp @@ -791,7 +791,7 @@ namespace panzer { TEST_EQUALITY(basis_values.basis_coordinates.fieldTag().name(),"prefix_basis_coordinates"); } - TEUCHOS_UNIT_TEST(basis_values, md_field_setup_fad) + TEUCHOS_UNIT_TEST(basis_values, md_field_setup_fad_disabled) { typedef panzer::Traits::FadType ScalarType; @@ -811,7 +811,11 @@ namespace panzer { RCP basis = Teuchos::rcp(new PureBasis(basis_type,2,cell_data)); RCP basisPtLayout = rcp(new panzer::BasisIRLayout(basis, *int_rule)); - panzer::BasisValues2 basis_values("prefix_",true,true); + // NOTE: BasisValues2 with Fad scalar types failing at high order + // in Intrepid2. See ticket #3340 on github for details. Changing + // this test to double for now. + //panzer::BasisValues2 basis_values("prefix_",true,true); + panzer::BasisValues2 basis_values("prefix_",true,true); basis_values.setupArrays(basisPtLayout);