diff --git a/packages/panzer/adapters-stk/example/PoissonInterfaceExample/main.cpp b/packages/panzer/adapters-stk/example/PoissonInterfaceExample/main.cpp index dd29f68ea5c5..822caf5ce5d2 100644 --- a/packages/panzer/adapters-stk/example/PoissonInterfaceExample/main.cpp +++ b/packages/panzer/adapters-stk/example/PoissonInterfaceExample/main.cpp @@ -62,6 +62,7 @@ #include "Panzer_BlockedEpetraLinearObjFactory.hpp" #include "Panzer_DOFManagerFactory.hpp" #include "Panzer_DOFManager.hpp" +#include "Panzer_OrientationsInterface.hpp" #include "Panzer_FieldManagerBuilder.hpp" #include "Panzer_PureBasis.hpp" #include "Panzer_GlobalData.hpp" @@ -568,23 +569,6 @@ int main (int argc, char* argv[]) mesh_factory->completeMeshConstruction(*mesh,MPI_COMM_WORLD); } - // build worksets - //////////////////////////////////////////////////////// - - Teuchos::RCP wkstFactory - = Teuchos::rcp(new panzer_stk::WorksetFactory(mesh)); // build STK workset factory - Teuchos::RCP wkstContainer // attach it to a workset container (uses lazy evaluation) - = Teuchos::rcp(new panzer::WorksetContainer); - wkstContainer->setFactory(wkstFactory); - for(size_t i=0;isetNeeds(physicsBlocks[i]->elementBlockID(),physicsBlocks[i]->getWorksetNeeds()); - wkstContainer->setWorksetSize(workset_size); - - std::vector elementBlockNames; - mesh->getElementBlockNames(elementBlockNames); - std::map > > volume_worksets; - panzer::getVolumeWorksetsFromContainer(*wkstContainer,elementBlockNames,volume_worksets); - // build DOF Manager and linear object factory ///////////////////////////////////////////////////////////// @@ -603,6 +587,25 @@ int main (int argc, char* argv[]) checkInterfaceConnections(conn_manager, dofManager->getComm()); } + // build worksets + //////////////////////////////////////////////////////// + + Teuchos::RCP wkstFactory + = Teuchos::rcp(new panzer_stk::WorksetFactory(mesh)); // build STK workset factory + wkstFactory->setOrientationsInterface(Teuchos::rcp(new panzer::OrientationsInterface(dofManager))); + Teuchos::RCP wkstContainer // attach it to a workset container (uses lazy evaluation) + = Teuchos::rcp(new panzer::WorksetContainer); + wkstContainer->setFactory(wkstFactory); + for(size_t i=0;isetNeeds(physicsBlocks[i]->elementBlockID(),physicsBlocks[i]->getWorksetNeeds()); + wkstContainer->setWorksetSize(workset_size); + + std::vector elementBlockNames; + mesh->getElementBlockNames(elementBlockNames); + std::map > > volume_worksets; + panzer::getVolumeWorksetsFromContainer(*wkstContainer,elementBlockNames,volume_worksets); + + // construct some linear algebra object, build object to pass to evaluators Teuchos::RCP > linObjFactory = Teuchos::rcp(new panzer::BlockedEpetraLinearObjFactory(tComm.getConst(),dofManager)); diff --git a/packages/panzer/adapters-stk/test/panzer_workset_builder/check_rotation_matrices.cpp b/packages/panzer/adapters-stk/test/panzer_workset_builder/check_rotation_matrices.cpp index 8e5e4f4271f8..4492c61b8039 100644 --- a/packages/panzer/adapters-stk/test/panzer_workset_builder/check_rotation_matrices.cpp +++ b/packages/panzer/adapters-stk/test/panzer_workset_builder/check_rotation_matrices.cpp @@ -120,13 +120,11 @@ namespace panzer { ////////////////////////////////////////////////////////////// panzer::IntegrationDescriptor sid(2*2, panzer::IntegrationDescriptor::SURFACE); - std::map wkstRequirements; - wkstRequirements[element_block].addIntegrator(sid); RCP wkstFactory = rcp(new panzer_stk::WorksetFactory(mesh)); // build STK workset factory RCP wkstContainer // attach it to a workset container (uses lazy evaluation) - = rcp(new panzer::WorksetContainer(wkstFactory,wkstRequirements)); + = rcp(new panzer::WorksetContainer(wkstFactory)); wkstContainer->setGlobalIndexer(dof_manager); @@ -138,8 +136,8 @@ namespace panzer { auto & workset = (*worksets)[0]; - auto rot_matrices = workset.getIntegrationValues(sid).surface_rotation_matrices; - auto normals = workset.getIntegrationValues(sid).surface_normals; + auto rot_matrices = workset.getIntegrationValues(sid).getSurfaceRotationMatrices(); + auto normals = workset.getIntegrationValues(sid).getSurfaceNormals(); const int num_owned_cells = workset.numOwnedCells(); const int num_ghost_cells = workset.numGhostCells(); diff --git a/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp b/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp index 20bbdc3e1fec..e34f4bd636e4 100644 --- a/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp +++ b/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp @@ -271,18 +271,18 @@ namespace panzer { bool is_uniform_; // Only valid for uniform reference space - PHX::MDField cubature_points_uniform_ref_; + PHX::MDField cubature_points_uniform_ref_; // For non-uniform reference space - PHX::MDField cubature_points_ref_; + PHX::MDField cubature_points_ref_; // Geometry objects that represent cubature/point values - PHX::MDField cubature_jacobian_; - PHX::MDField cubature_jacobian_determinant_; - PHX::MDField cubature_jacobian_inverse_; - PHX::MDField cubature_weights_; + PHX::MDField cubature_jacobian_; + PHX::MDField cubature_jacobian_determinant_; + PHX::MDField cubature_jacobian_inverse_; + PHX::MDField cubature_weights_; - PHX::MDField cell_vertex_coordinates_; + PHX::MDField cell_vertex_coordinates_; // Number of cells to apply orientations to (required in situations where virtual cells exist) int num_orientations_cells_; @@ -329,10 +329,10 @@ namespace panzer { */ void setup(const Teuchos::RCP & basis, - PHX::MDField reference_points, - PHX::MDField point_jacobian, - PHX::MDField point_jacobian_determinant, - PHX::MDField point_jacobian_inverse, + PHX::MDField reference_points, + PHX::MDField point_jacobian, + PHX::MDField point_jacobian_determinant, + PHX::MDField point_jacobian_inverse, const int num_evaluated_cells = -1); /** @@ -350,10 +350,10 @@ namespace panzer { */ void setupUniform(const Teuchos::RCP & basis, - PHX::MDField reference_points, - PHX::MDField point_jacobian, - PHX::MDField point_jacobian_determinant, - PHX::MDField point_jacobian_inverse, + PHX::MDField reference_points, + PHX::MDField point_jacobian, + PHX::MDField point_jacobian_determinant, + PHX::MDField point_jacobian_inverse, const int num_evaluated_cells = -1); /// Set the orientations object for applying orientations using the lazy evaluation path - required for certain bases @@ -363,7 +363,7 @@ namespace panzer { /// Set the cubature weights (weighted measure) for the basis values object - required to get weighted basis objects void - setWeightedMeasure(PHX::MDField weighted_measure); + setWeightedMeasure(PHX::MDField weighted_measure); /// Set the cell vertex coordinates (required for getBasisCoordinates()) void diff --git a/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp b/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp index fcfd580eb29b..c4fdecebc00b 100644 --- a/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp +++ b/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp @@ -52,6 +52,9 @@ #include "Intrepid2_Orientation.hpp" #include "Intrepid2_OrientationTools.hpp" +// FIXME: There are some calls in Intrepid2 that require non-const arrays when they should be const - search for PHX::getNonConstDynRankViewFromConstMDField +#include "Phalanx_GetNonConstDynRankViewFromConstMDField.hpp" + namespace panzer { namespace { @@ -599,10 +602,10 @@ template void BasisValues2:: setup(const Teuchos::RCP & basis, - PHX::MDField reference_points, - PHX::MDField point_jacobian, - PHX::MDField point_jacobian_determinant, - PHX::MDField point_jacobian_inverse, + PHX::MDField reference_points, + PHX::MDField point_jacobian, + PHX::MDField point_jacobian_determinant, + PHX::MDField point_jacobian_inverse, const int num_evaluated_cells) { basis_layout = basis; @@ -626,10 +629,10 @@ template void BasisValues2:: setupUniform(const Teuchos::RCP & basis, - PHX::MDField reference_points, - PHX::MDField point_jacobian, - PHX::MDField point_jacobian_determinant, - PHX::MDField point_jacobian_inverse, + PHX::MDField reference_points, + PHX::MDField point_jacobian, + PHX::MDField point_jacobian_determinant, + PHX::MDField point_jacobian_inverse, const int num_evaluated_cells) { basis_layout = basis; @@ -735,7 +738,7 @@ resetArrays() template void BasisValues2:: -setWeightedMeasure(PHX::MDField weighted_measure) +setWeightedMeasure(PHX::MDField weighted_measure) { TEUCHOS_TEST_FOR_EXCEPT_MSG(build_weighted, "BasisValues2::setWeightedMeasure : Weighted measure already set. Can only set weighted measure once after setup or setupUniform have beens called."); @@ -807,7 +810,9 @@ getBasisValuesRef(const bool cache, const int num_card = basis_layout->cardinality(); auto tmp_basis_ref_scalar = af.buildStaticArray("dyn_basis_ref_scalar",num_card,num_quad); - intrepid_basis->getValues(tmp_basis_ref_scalar.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_VALUE); + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + + intrepid_basis->getValues(tmp_basis_ref_scalar.get_view(), cubature_points_uniform_ref, Intrepid2::OPERATOR_VALUE); PHX::Device().fence(); // Store for later if cache is enabled @@ -840,7 +845,9 @@ getVectorBasisValuesRef(const bool cache, const int num_dim = basis_layout->dimension(); auto tmp_basis_ref_vector = af.buildStaticArray("dyn_basis_ref_vector",num_card,num_quad,num_dim); - intrepid_basis->getValues(tmp_basis_ref_vector.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_VALUE); + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + + intrepid_basis->getValues(tmp_basis_ref_vector.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_VALUE); PHX::Device().fence(); // Store for later if cache is enabled @@ -873,7 +880,9 @@ getGradBasisValuesRef(const bool cache, const int num_dim = basis_layout->dimension(); auto tmp_grad_basis_ref = af.buildStaticArray("dyn_basis_ref_vector",num_card,num_quad,num_dim); - intrepid_basis->getValues(tmp_grad_basis_ref.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_GRAD); + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + + intrepid_basis->getValues(tmp_grad_basis_ref.get_view(), cubature_points_uniform_ref, Intrepid2::OPERATOR_GRAD); PHX::Device().fence(); // Store for later if cache is enabled @@ -906,7 +915,9 @@ getCurl2DVectorBasisRef(const bool cache, const int num_card = basis_layout->cardinality(); auto tmp_curl_basis_ref_scalar = af.buildStaticArray("dyn_curl_basis_ref_scalar",num_card,num_quad); - intrepid_basis->getValues(tmp_curl_basis_ref_scalar.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_CURL); + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + + intrepid_basis->getValues(tmp_curl_basis_ref_scalar.get_view(), cubature_points_uniform_ref, Intrepid2::OPERATOR_CURL); PHX::Device().fence(); // Store for later if cache is enabled @@ -940,7 +951,9 @@ getCurlVectorBasisRef(const bool cache, const int num_dim = basis_layout->dimension(); auto tmp_curl_basis_ref_vector = af.buildStaticArray("dyn_curl_basis_ref_vector",num_card,num_quad,num_dim); - intrepid_basis->getValues(tmp_curl_basis_ref_vector.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_CURL); + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + + intrepid_basis->getValues(tmp_curl_basis_ref_vector.get_view(), cubature_points_uniform_ref, Intrepid2::OPERATOR_CURL); PHX::Device().fence(); // Store for later if cache is enabled @@ -972,7 +985,9 @@ getDivVectorBasisRef(const bool cache, const int num_card = basis_layout->cardinality(); auto tmp_div_basis_ref = af.buildStaticArray("dyn_div_basis_ref_scalar",num_card,num_quad); - intrepid_basis->getValues(tmp_div_basis_ref.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_DIV); + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + + intrepid_basis->getValues(tmp_div_basis_ref.get_view(), cubature_points_uniform_ref, Intrepid2::OPERATOR_DIV); PHX::Device().fence(); // Store for later if cache is enabled @@ -1004,13 +1019,12 @@ getBasisCoordinates(const bool cache, auto tmp_basis_coordinates = af.buildStaticArray("basis_coordinates", num_cells_, num_card, num_dim); auto s_aux = Kokkos::subview(tmp_basis_coordinates.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL()); - // Some kind of type bug in Intrepid's mapToPhysicalFrame call? Can't let bcr be const - auto bcr = getBasisCoordinatesRef(false); - Kokkos::DynRankView bcr_copy("nonconst_bcr",bcr.extent(0),bcr.extent(1)); - Kokkos::deep_copy(bcr_copy,bcr.get_view()); + // Don't forget that since we are not caching this, we have to make sure the managed view remains alive while we use the non-const wrapper + auto const_bcr = getBasisCoordinatesRef(false); + auto bcr = PHX::getNonConstDynRankViewFromConstMDField(const_bcr); Intrepid2::CellTools cell_tools; - cell_tools.mapToPhysicalFrame(s_aux, bcr_copy, s_vertex_coordinates, intrepid_basis->getBaseCellTopology()); + cell_tools.mapToPhysicalFrame(s_aux, bcr, s_vertex_coordinates, intrepid_basis->getBaseCellTopology()); PHX::Device().fence(); // Store for later if cache is enabled @@ -1081,8 +1095,10 @@ getBasisValues(const bool weighted, if(hasUniformReferenceSpace()){ + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + // Apply a single reference representation to all cells - intrepid_basis->getValues(cell_basis_ref_scalar.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_VALUE); + intrepid_basis->getValues(cell_basis_ref_scalar.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_VALUE); const std::pair cell_range(0,num_evaluate_cells_); auto s_aux = Kokkos::subview(tmp_basis_scalar.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL()); @@ -1239,8 +1255,10 @@ getVectorBasisValues(const bool weighted, if(hasUniformReferenceSpace()){ + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + // Apply a single reference representation to all cells - intrepid_basis->getValues(cell_basis_ref_vector.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_VALUE); + intrepid_basis->getValues(cell_basis_ref_vector.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_VALUE); const std::pair cell_range(0,num_evaluate_cells_); auto s_aux = Kokkos::subview(tmp_basis_vector.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); @@ -1414,8 +1432,10 @@ getGradBasisValues(const bool weighted, if(hasUniformReferenceSpace()){ + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + // Apply a single reference representation to all cells - intrepid_basis->getValues(cell_grad_basis_ref.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_GRAD); + intrepid_basis->getValues(cell_grad_basis_ref.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_GRAD); const std::pair cell_range(0,num_evaluate_cells_); auto s_aux = Kokkos::subview(tmp_grad_basis.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); @@ -1560,7 +1580,9 @@ getCurl2DVectorBasis(const bool weighted, if(hasUniformReferenceSpace()){ - intrepid_basis->getValues(cell_curl_basis_ref_scalar.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_CURL); + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + + intrepid_basis->getValues(cell_curl_basis_ref_scalar.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_CURL); const std::pair cell_range(0,num_evaluate_cells_); auto s_aux = Kokkos::subview(tmp_curl_basis_scalar.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL()); @@ -1703,7 +1725,9 @@ getCurlVectorBasis(const bool weighted, if(hasUniformReferenceSpace()){ - intrepid_basis->getValues(cell_curl_basis_ref_vector.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_CURL); + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + + intrepid_basis->getValues(cell_curl_basis_ref_vector.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_CURL); const std::pair cell_range(0,num_evaluate_cells_); auto s_aux = Kokkos::subview(tmp_curl_basis_vector.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); @@ -1848,7 +1872,9 @@ getDivVectorBasis(const bool weighted, if(hasUniformReferenceSpace()){ - intrepid_basis->getValues(cell_div_basis_ref.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_DIV); + auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_); + + intrepid_basis->getValues(cell_div_basis_ref.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_DIV); const std::pair cell_range(0,num_evaluate_cells_); auto s_aux = Kokkos::subview(tmp_div_basis.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL()); diff --git a/packages/panzer/disc-fe/src/Panzer_IntegrationDescriptor.cpp b/packages/panzer/disc-fe/src/Panzer_IntegrationDescriptor.cpp index 9732e3d39aa0..422160b775ab 100644 --- a/packages/panzer/disc-fe/src/Panzer_IntegrationDescriptor.cpp +++ b/packages/panzer/disc-fe/src/Panzer_IntegrationDescriptor.cpp @@ -66,7 +66,7 @@ IntegrationDescriptor::setup(const int cubature_order, const int integration_typ _cubature_order = cubature_order; _side = side; - if(_integration_type == SIDE or _integration_type == CV_BOUNDARY){ + if((_integration_type == SIDE) or (_integration_type == CV_BOUNDARY)){ TEUCHOS_ASSERT(side >= 0); } else { TEUCHOS_ASSERT(side == -1); diff --git a/packages/panzer/disc-fe/src/Panzer_IntegrationRule.cpp b/packages/panzer/disc-fe/src/Panzer_IntegrationRule.cpp index 475504328806..a07b6d39a792 100644 --- a/packages/panzer/disc-fe/src/Panzer_IntegrationRule.cpp +++ b/packages/panzer/disc-fe/src/Panzer_IntegrationRule.cpp @@ -72,25 +72,16 @@ IntegrationRule(const panzer::CellData& cell_data, const std::string & in_cv_typ { // Cubature orders are only used for indexing so we make them large enough not to interfere with other rules. // TODO: This requirement (on arbitrary cubature order) will be dropped with the new Workset design (using descriptors to find integration rules) - // TODO: These isSide conditions shouldn't be required... I'm missing something... if(in_cv_type == "volume"){ - if(cell_data.isSide()){ - IntegrationDescriptor::setup(75, IntegrationDescriptor::CV_VOLUME,cell_data.side()); - } else { - IntegrationDescriptor::setup(75, IntegrationDescriptor::CV_VOLUME); - } + TEUCHOS_TEST_FOR_EXCEPT_MSG(cell_data.isSide(), + "IntegrationRule::IntegrationRule : Control Volume 'volume' type requested, but CellData is setup for sides."); + IntegrationDescriptor::setup(75, IntegrationDescriptor::CV_VOLUME); } else if(in_cv_type == "side"){ - if(cell_data.isSide()){ - IntegrationDescriptor::setup(85, IntegrationDescriptor::CV_SIDE,cell_data.side()); - } else { - IntegrationDescriptor::setup(85, IntegrationDescriptor::CV_SIDE); - } + IntegrationDescriptor::setup(85, IntegrationDescriptor::CV_SIDE); } else if(in_cv_type == "boundary"){ - if(cell_data.isSide()){ - IntegrationDescriptor::setup(95, IntegrationDescriptor::CV_BOUNDARY, cell_data.side()); - } else { - IntegrationDescriptor::setup(95, IntegrationDescriptor::CV_BOUNDARY); - } + TEUCHOS_TEST_FOR_EXCEPT_MSG(not cell_data.isSide(), + "IntegrationRule::IntegrationRule : Control Volume 'boundary' type requested, but CellData is not setup for sides."); + IntegrationDescriptor::setup(95, IntegrationDescriptor::CV_BOUNDARY, cell_data.side()); } else { TEUCHOS_ASSERT(false); } diff --git a/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.cpp b/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.cpp index f8a5b9e807b4..cbe5c45a2d3d 100644 --- a/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.cpp +++ b/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.cpp @@ -59,179 +59,15 @@ #include "Panzer_SubcellConnectivity.hpp" #include "Panzer_ConvertNormalToRotationMatrix.hpp" -// For device member functions to avoid CUDA RDC in unit testing -#include "Panzer_IntegrationValues2_impl.hpp" - -// *********************************************************** -// * Specializations of setupArrays() for different array types -// *********************************************************** +// FIXME: There are some calls in Intrepid2 that require non-const arrays when they should be const - search for PHX::getNonConstDynRankViewFromConstMDField +#include "Phalanx_GetNonConstDynRankViewFromConstMDField.hpp" namespace panzer { -// * Specialization for Kokkos::DynRankView -template -void IntegrationValues2:: -setupArraysForNodeRule(const Teuchos::RCP& ir) -{ - MDFieldArrayFactory af(prefix,alloc_arrays); - - int num_nodes = ir->topology->getNodeCount(); - int num_cells = ir->workset_size; - int num_space_dim = ir->topology->getDimension(); - - int num_ip = 1; - - dyn_cub_points = af.template buildArray("cub_points",num_ip, num_space_dim); - dyn_cub_weights = af.template buildArray("cub_weights",num_ip); - - cub_points = af.template buildStaticArray("cub_points",num_ip, num_space_dim); - - if (ir->cv_type == "none" && ir->isSide()) { - dyn_side_cub_points = af.template buildArray("side_cub_points",num_ip, ir->side_topology->getDimension()); - side_cub_points = af.template buildStaticArray("side_cub_points",num_ip,ir->side_topology->getDimension()); - } - - if (ir->cv_type != "none") { - dyn_phys_cub_points = af.template buildArray("phys_cub_points",num_cells, num_ip, num_space_dim); - dyn_phys_cub_weights = af.template buildArray("phys_cub_weights",num_cells, num_ip); - if (ir->cv_type == "side") { - dyn_phys_cub_norms = af.template buildArray("phys_cub_norms",num_cells, num_ip, num_space_dim); - } - } - - dyn_node_coordinates = af.template buildArray("node_coordinates",num_cells, num_nodes, num_space_dim); - - cub_weights = af.template buildStaticArray("cub_weights",num_ip); - - node_coordinates = af.template buildStaticArray("node_coordinates",num_cells, num_nodes, num_space_dim); - - jac = af.template buildStaticArray("jac",num_cells, num_ip, num_space_dim,num_space_dim); - - jac_inv = af.template buildStaticArray("jac_inv",num_cells, num_ip, num_space_dim,num_space_dim); - - jac_det = af.template buildStaticArray("jac_det",num_cells, num_ip); - - weighted_measure = af.template buildStaticArray("weighted_measure",num_cells, num_ip); - - covarient = af.template buildStaticArray("covarient",num_cells, num_ip, num_space_dim,num_space_dim); - - contravarient = af.template buildStaticArray("contravarient",num_cells, num_ip, num_space_dim,num_space_dim); - - norm_contravarient = af.template buildStaticArray("norm_contravarient",num_cells, num_ip); - - ip_coordinates = af.template buildStaticArray("ip_coordiantes",num_cells, num_ip,num_space_dim); - - ref_ip_coordinates = af.template buildStaticArray("ref_ip_coordinates",num_cells, num_ip,num_space_dim); - - weighted_normals = af.template buildStaticArray("weighted normal",num_cells, num_ip,num_space_dim); - - surface_normals = af.template buildStaticArray("surface_normals",num_cells, num_ip,num_space_dim); - - surface_rotation_matrices = af.template buildStaticArray("surface_rotation_matrices",num_cells, num_ip,3,3); - -} - -template -void IntegrationValues2:: -setupArrays(const Teuchos::RCP& ir) -{ - MDFieldArrayFactory af(prefix,alloc_arrays); - - typedef panzer::IntegrationDescriptor ID; - - int_rule = ir; - - int num_nodes = ir->topology->getNodeCount(); - int num_cells = ir->workset_size; - int num_space_dim = ir->topology->getDimension(); - - // specialize content if this is quadrature at anode - if(num_space_dim==1 && ir->isSide()) { - setupArraysForNodeRule(ir); - return; - } - - TEUCHOS_ASSERT(ir->getType() != ID::NONE); - intrepid_cubature = getIntrepidCubature(*ir); - - int num_ip = ir->num_points; - -// Intrepid2::DefaultCubatureFactory cubature_factory; -// -// if (ir->cv_type == "side") -// intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolumeSide(*ir->topology)); -// -// else if (ir->cv_type == "volume") -// intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolume(*ir->topology)); -// -// else if (ir->cv_type == "boundary" && ir->isSide()) -// intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolumeBoundary(*ir->topology,ir->side)); -// -// else if (ir->cv_type == "none" && ir->isSide()) -// intrepid_cubature = cubature_factory.create(*(ir->side_topology), -// ir->cubature_degree); -// else -// intrepid_cubature = cubature_factory.create(*(ir->topology), -// ir->cubature_degree); -// int num_ip = intrepid_cubature->getNumPoints(); - - dyn_cub_points = af.template buildArray("cub_points",num_ip, num_space_dim); - dyn_cub_weights = af.template buildArray("cub_weights",num_ip); - - cub_points = af.template buildStaticArray("cub_points",num_ip, num_space_dim); - - if (ir->isSide() && ir->cv_type == "none") { - dyn_side_cub_points = af.template buildArray("side_cub_points",num_ip, ir->side_topology->getDimension()); - side_cub_points = af.template buildStaticArray("side_cub_points",num_ip,ir->side_topology->getDimension()); - } - - if (ir->cv_type != "none") { - dyn_phys_cub_points = af.template buildArray("phys_cub_points",num_cells, num_ip, num_space_dim); - dyn_phys_cub_weights = af.template buildArray("phys_cub_weights",num_cells, num_ip); - if (ir->cv_type == "side") { - dyn_phys_cub_norms = af.template buildArray("phys_cub_norms",num_cells, num_ip, num_space_dim); - } - } - - dyn_node_coordinates = af.template buildArray("node_coordinates",num_cells, num_nodes, num_space_dim); - - cub_weights = af.template buildStaticArray("cub_weights",num_ip); - - node_coordinates = af.template buildStaticArray("node_coordinates",num_cells, num_nodes, num_space_dim); - - jac = af.template buildStaticArray("jac",num_cells, num_ip, num_space_dim,num_space_dim); - - jac_inv = af.template buildStaticArray("jac_inv",num_cells, num_ip, num_space_dim,num_space_dim); - - jac_det = af.template buildStaticArray("jac_det",num_cells, num_ip); - - weighted_measure = af.template buildStaticArray("weighted_measure",num_cells, num_ip); - - covarient = af.template buildStaticArray("covarient",num_cells, num_ip, num_space_dim,num_space_dim); - - contravarient = af.template buildStaticArray("contravarient",num_cells, num_ip, num_space_dim,num_space_dim); - - norm_contravarient = af.template buildStaticArray("norm_contravarient",num_cells, num_ip); - - ip_coordinates = af.template buildStaticArray("ip_coordiantes",num_cells, num_ip,num_space_dim); - - ref_ip_coordinates = af.template buildStaticArray("ref_ip_coordinates",num_cells, num_ip,num_space_dim); - - weighted_normals = af.template buildStaticArray("weighted_normal",num_cells,num_ip,num_space_dim); - - surface_normals = af.template buildStaticArray("surface_normals",num_cells, num_ip,num_space_dim); - - surface_rotation_matrices = af.template buildStaticArray("surface_rotation_matrices",num_cells, num_ip,3,3); - - scratch_for_compute_side_measure = - af.template buildStaticArray("scratch_for_compute_side_measure", jac.get_view().span()); - -} +namespace { -template Teuchos::RCP> -IntegrationValues2:: -getIntrepidCubature(const panzer::IntegrationRule & ir) const +getIntrepidCubature(const panzer::IntegrationRule & ir) { typedef panzer::IntegrationDescriptor ID; Teuchos::RCP > ic; @@ -243,335 +79,333 @@ getIntrepidCubature(const panzer::IntegrationRule & ir) const } else if(ir.getType() == ID::CV_VOLUME){ ic = Teuchos::rcp(new Intrepid2::CubatureControlVolume(*ir.topology)); } else if(ir.getType() == ID::CV_BOUNDARY){ + TEUCHOS_ASSERT(ir.isSide()); ic = Teuchos::rcp(new Intrepid2::CubatureControlVolumeBoundary(*ir.topology,ir.getSide())); + } else if(ir.getType() == ID::VOLUME){ + ic = cubature_factory.create(*(ir.topology),ir.getOrder()); + } else if(ir.getType() == ID::SIDE){ + ic = cubature_factory.create(*(ir.side_topology),ir.getOrder()); + } else if(ir.getType() == ID::SURFACE){ + // closed surface integrals don't exist in intrepid. } else { - if(ir.getType() == ID::VOLUME){ - ic = cubature_factory.create(*(ir.topology),ir.getOrder()); - } else if(ir.getType() == ID::SIDE){ - ic = cubature_factory.create(*(ir.side_topology),ir.getOrder()); - } else if(ir.getType() == ID::SURFACE){ - // closed surface integrals don't exist in intrepid. - } else { - TEUCHOS_ASSERT(false); - } + TEUCHOS_ASSERT(false); } return ic; } - -// *********************************************************** -// * Evaluation of values - NOT specialized -// *********************************************************** -template -void IntegrationValues2:: -evaluateValues(const PHX::MDField & in_node_coordinates, - const int in_num_cells, - const Teuchos::RCP & face_connectivity) -{ - typedef panzer::IntegrationDescriptor ID; - const bool is_surface = int_rule->getType() == ID::SURFACE; - const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY); - - TEUCHOS_ASSERT(not (is_surface and is_cv)); - - if(is_surface){ - TEUCHOS_TEST_FOR_EXCEPT_MSG(face_connectivity == Teuchos::null, - "IntegrationValues2::evaluateValues : Surface integration requires the face connectivity"); - generateSurfaceCubatureValues(in_node_coordinates,in_num_cells,*face_connectivity); - } else if (is_cv) { - getCubatureCV(in_node_coordinates, in_num_cells); - evaluateValuesCV(in_node_coordinates, in_num_cells); - } else { - getCubature(in_node_coordinates, in_num_cells); - evaluateRemainingValues(in_node_coordinates, in_num_cells); - } -} - -template -void IntegrationValues2:: -getCubature(const PHX::MDField& in_node_coordinates, - const int in_num_cells) -{ - - int num_space_dim = int_rule->topology->getDimension(); - if (int_rule->isSide() && num_space_dim==1) { - std::cout << "WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do " - << "non-natural integration rules."; - return; - } - - Intrepid2::CellTools cell_tools; - - if (!int_rule->isSide()) - intrepid_cubature->getCubature(dyn_cub_points.get_view(), dyn_cub_weights.get_view()); - else { - intrepid_cubature->getCubature(dyn_side_cub_points.get_view(), dyn_cub_weights.get_view()); - - cell_tools.mapToReferenceSubcell(dyn_cub_points.get_view(), - dyn_side_cub_points.get_view(), - int_rule->spatial_dimension-1, - int_rule->side, - *(int_rule->topology)); - } - - // IP coordinates - const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells; - auto s_ip_coordinates = Kokkos::subview(ip_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL()); - auto s_in_node_coordinates = Kokkos::subview(in_node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL()); - cell_tools.mapToPhysicalFrame(s_ip_coordinates, - dyn_cub_points.get_view(), - s_in_node_coordinates, - *(int_rule->topology)); -} - -template -void IntegrationValues2:: -generateSurfaceCubatureValues(const PHX::MDField& in_node_coordinates, - const int in_num_cells, - const SubcellConnectivity & face_connectivity) +template +void +correctVirtualNormals(PHX::MDField normals, + const int num_real_cells, + const int num_virtual_cells, + const shards::CellTopology & cell_topology, + const SubcellConnectivity & face_connectivity) { - TEUCHOS_ASSERT(int_rule->getType() == IntegrationDescriptor::SURFACE); - - Intrepid2::CellTools cell_tools; - - const shards::CellTopology & cell_topology = *(int_rule->topology); - const panzer::IntegrationRule & ir = *int_rule; - - const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells; - - // Copy over coordinates - { - const int num_nodes = in_node_coordinates.extent(1); - const int num_dims = in_node_coordinates.extent(2); - auto node_coordinates_k = node_coordinates.get_view(); - auto in_node_coordinates_k = in_node_coordinates.get_view(); - - Kokkos::MDRangePolicy> policy({0,0,0},{num_cells,num_nodes,num_dims}); - Kokkos::parallel_for("node_coordinates",policy,KOKKOS_LAMBDA (const int cell, const int node, const int dim) { - node_coordinates_k(cell,node,dim) = in_node_coordinates_k(cell,node,dim); - }); - PHX::Device::execution_space().fence(); - } + // What we want is for the adjoining face of the virtual cell to have normals that are the negated real cell's normals. + // we correct the normals here: + const int space_dim = cell_topology.getDimension(); + const int faces_per_cell = cell_topology.getSubcellCount(space_dim-1); + const int points = normals.extent_int(1); + const int points_per_face = points / faces_per_cell; + + Kokkos::parallel_for(num_virtual_cells, KOKKOS_LAMBDA(const int & virtual_cell_ordinal){ + + const int virtual_cell = virtual_cell_ordinal+num_real_cells; + + // Find the face and local face ids for the given virtual cell + // Note that virtual cells only connect to the owned cells through a single face + int face, virtual_lidx; + for (int local_face_id=0; local_face_id= 0){ + virtual_lidx = local_face_id; + break; + } + } - // NOTE: We are assuming that each face can have a different number of points. - // Not sure if this is necessary, but it requires a lot of additional allocations + // Indexes for real cell + const int real_side = (face_connectivity.cellForSubcell(face, 0) == virtual_cell) ? 1 : 0; + const int real_cell = face_connectivity.cellForSubcell(face,real_side); + const int real_lidx = face_connectivity.localSubcellForSubcell(face,real_side); - const int cell_dim = cell_topology.getDimension(); - const int subcell_dim = cell_topology.getDimension()-1; - const int num_subcells = cell_topology.getSubcellCount(subcell_dim); + // Make sure it is a real cell (it should actually be an owned cell) + KOKKOS_ASSERT(real_cell < num_real_cells); - Intrepid2::DefaultCubatureFactory cubature_factory; + for(int point_ordinal=0; point_ordinal tmp_side_cub_weights; - Kokkos::DynRankView tmp_side_cub_points; - if(cell_dim==1){ - tmp_side_cub_weights = Kokkos::DynRankView("tmp_side_cub_weights",num_points_on_face); - tmp_side_cub_points = Kokkos::DynRankView("cell_tmp_side_cub_points",num_points_on_face,cell_dim); - auto tmp_side_cub_weights_host = Kokkos::create_mirror_view(tmp_side_cub_weights); - auto tmp_side_cub_points_host = Kokkos::create_mirror_view(tmp_side_cub_points); - tmp_side_cub_weights_host(0)=1.; - tmp_side_cub_points_host(0,0) = (subcell_index==0)? -1. : 1.; - Kokkos::deep_copy(tmp_side_cub_weights,tmp_side_cub_weights_host); - Kokkos::deep_copy(tmp_side_cub_points,tmp_side_cub_points_host); - } else { + } - // Get the face topology from the cell topology - const shards::CellTopology face_topology(cell_topology.getCellTopologyData(subcell_dim,subcell_index)); + // Clear other normals + for (int local_face_id=0; local_face_id(face_topology,ir.getOrder()); - num_points_on_face = ic->getNumPoints(); + if (local_face_id == virtual_lidx) continue; - tmp_side_cub_weights = Kokkos::DynRankView("tmp_side_cub_weights",num_points_on_face); - tmp_side_cub_points = Kokkos::DynRankView("cell_tmp_side_cub_points",num_points_on_face,cell_dim); + for (int point_ordinal=0; point_ordinal("subcell_cub_points",num_points_on_face,subcell_dim); - // Get the reference face points - ic->getCubature(subcell_cub_points, tmp_side_cub_weights); - - // Convert from reference face points to reference cell points - cell_tools.mapToReferenceSubcell(tmp_side_cub_points, - subcell_cub_points, - subcell_dim, - subcell_index, - cell_topology); - } +template +void +correctVirtualRotationMatrices(PHX::MDField rotation_matrices, + const int num_real_cells, + const int num_virtual_cells, + const shards::CellTopology & cell_topology, + const SubcellConnectivity & face_connectivity) +{ - // Do this on host - { - auto tmp_side_cub_points_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),tmp_side_cub_points); - auto cub_points_host = Kokkos::create_mirror_view(cub_points.get_static_view()); - for(int local_point=0;local_point= 0){ + virtual_lidx = local_face_id; + break; } + } + // The normals already have the correction applied, so we just need to zero out the rotation matrices on the other faces - // Map from side points to physical points - auto side_ip_coordinates = Kokkos::DynRankView("side_ip_coordinates",num_cells,num_points_on_face,cell_dim); - auto s_node_coordinates = Kokkos::subview(node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL); - cell_tools.mapToPhysicalFrame(side_ip_coordinates, - tmp_side_cub_points, - s_node_coordinates, - cell_topology); - - // Create a jacobian and his friends for this side - auto side_jacobian = Kokkos::DynRankView("side_jac",num_cells,num_points_on_face,cell_dim,cell_dim); - cell_tools.setJacobian(side_jacobian, - tmp_side_cub_points, - s_node_coordinates, - cell_topology); + // Clear other rotation matrices + for (int local_face_id=0; local_face_id("side_inv_jac",num_cells,num_points_on_face,cell_dim,cell_dim); - cell_tools.setJacobianInv(side_inverse_jacobian, side_jacobian); + for (int point_ordinal=0; point_ordinal("side_det_jac",num_cells,num_points_on_face); - cell_tools.setJacobianDet(side_det_jacobian, side_jacobian); +template +void +applyBasePermutation(PHX::MDField field, + PHX::MDField permutations) +{ + MDFieldArrayFactory af("",true); - PHX::Device::execution_space().fence(); + const int num_ip = field.extent(0); - // Calculate measures (quadrature weights in physical space) for this side - auto side_weighted_measure = Kokkos::DynRankView("side_weighted_measure",num_cells,num_points_on_face); - if(cell_dim == 1){ - Kokkos::deep_copy(side_weighted_measure, tmp_side_cub_weights(0)); - } else if(cell_dim == 2){ - Intrepid2::FunctionSpaceTools:: - computeEdgeMeasure(side_weighted_measure, side_jacobian, tmp_side_cub_weights, - subcell_index,cell_topology, - scratch_for_compute_side_measure.get_view()); - PHX::Device::execution_space().fence(); - } else if(cell_dim == 3){ - Intrepid2::FunctionSpaceTools:: - computeFaceMeasure(side_weighted_measure, side_jacobian, tmp_side_cub_weights, - subcell_index,cell_topology, - scratch_for_compute_side_measure.get_view()); - PHX::Device::execution_space().fence(); - } + auto scratch = af.template buildStaticArray("scratch", num_ip); + scratch.deep_copy(field); + Kokkos::parallel_for(1, KOKKOS_LAMBDA(const int & ){ + for(int ip=0; ip("side_normals",num_cells,num_points_on_face,cell_dim); - if(cell_dim == 1){ +template +void +applyBasePermutation(PHX::MDField field, + PHX::MDField permutations) +{ + MDFieldArrayFactory af("",true); + + const int num_ip = field.extent(0); + const int num_dim = field.extent(1); + + auto scratch = af.template buildStaticArray("scratch", num_ip,num_dim); + scratch.deep_copy(field); + Kokkos::parallel_for(1, KOKKOS_LAMBDA(const int & ){ + for(int ip=0; ip::type>::min(); +template +void +applyPermutation(PHX::MDField field, + PHX::MDField permutations) +{ + MDFieldArrayFactory af("",true); - Kokkos::parallel_for("compute normals 1D",num_cells,KOKKOS_LAMBDA (const int cell) { - Scalar norm = (in_node_coordinates_k(cell,subcell_index,0) - in_node_coordinates_k(cell,other_subcell_index,0)); - side_normals(cell,0,0) = norm / fabs(norm + min); - }); + const int num_cells = field.extent(0); + const int num_ip = field.extent(1); - } else { + auto scratch = af.template buildStaticArray("scratch", num_cells, num_ip); + scratch.deep_copy(field); + Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(const int & cell){ + for(int ip=0; ip +void +applyPermutation(PHX::MDField field, + PHX::MDField permutations) +{ + MDFieldArrayFactory af("",true); + + const int num_cells = field.extent(0); + const int num_ip = field.extent(1); + const int num_dim = field.extent(2); + + auto scratch = af.template buildStaticArray("scratch", num_cells, num_ip, num_dim); + scratch.deep_copy(field); + Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(const int & cell){ + for(int ip=0; ip> policy({0,0},{num_cells,num_points_on_face}); - Kokkos::parallel_for("Normalize the normals",policy,KOKKOS_LAMBDA (const int cell,const int point) { - Scalar n = 0.; - for(int dim=0;dim 0.){ - n = Kokkos::Experimental::sqrt(n); - for(int dim=0;dim +void +applyPermutation(PHX::MDField field, + PHX::MDField permutations) +{ + MDFieldArrayFactory af("",true); + + const int num_cells = field.extent(0); + const int num_ip = field.extent(1); + const int num_dim = field.extent(2); + const int num_dim2 = field.extent(3); + + auto scratch = af.template buildStaticArray("scratch", num_cells, num_ip, num_dim, num_dim2); + scratch.deep_copy(field); + Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(const int & cell){ + for(int ip=0; ip> policy({0,0},{num_cells,num_points_on_face}); - Kokkos::parallel_for("copy values",policy,KOKKOS_LAMBDA (const int cell,const int side_point) { - const int cell_point = point_offset + side_point; - - weighted_measure_k(cell,cell_point) = side_weighted_measure(cell,side_point); - jac_det_k(cell,cell_point) = side_det_jacobian(cell,side_point); - for(int dim=0;dim +PHX::MDField +generatePermutations(const int num_cells, + PHX::MDField coords, + PHX::MDField other_coords) +{ + const int num_ip = coords.extent(1); + const int num_dim = coords.extent(2); + + MDFieldArrayFactory af("",true); + + // This will store the permutations to go from coords to other_coords + auto permutation = af.template buildStaticArray("permutation", num_cells, num_ip); + permutation.deep_copy(0); + + // This is scratch space - there is likely a better way to do this + auto taken = af.template buildStaticArray("taken", num_cells, num_ip); + taken.deep_copy(0); + + Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(const int & cell){ + + for (int ip = 0; ip < num_ip; ++ip) { + // Find an other point to associate with ip. + int i_min = 0; + Scalar d_min = -1; + for (int other_ip = 0; other_ip < num_ip; ++other_ip) { + // For speed, skip other points that are already associated. + if(taken(cell,other_ip)) continue; + // Compute the distance between the two points. + Scalar d(0); + for (int dim = 0; dim < num_dim; ++dim) { + const Scalar diff = coords(cell, ip, dim) - other_coords(cell, other_ip, dim); + d += diff*diff; + } + if (d_min < 0 || d < d_min) { + d_min = d; + i_min = other_ip; + } } - point_offset += num_points_on_face; + // Record the match. + permutation(cell,ip) = i_min; + // This point is now taken. + taken(cell,i_min) = 1; } - } - - // We also need surface rotation matrices in order to enforce alignment - { - const int num_points = ir.getPointOffset(num_subcells); - auto surface_normals_k = surface_normals.get_view(); - auto surface_rotation_matrices_k = surface_rotation_matrices.get_view(); - Kokkos::MDRangePolicy> policy({0,0,0},{num_cells,num_subcells,num_points}); - Kokkos::parallel_for("create surface rotation matrices",policy,KOKKOS_LAMBDA (const int cell,const int subcell_index,const int point) { - Scalar normal[3]; + }); + PHX::Device::execution_space().fence(); - for(int i=0;i<3;i++){normal[i]=0.;} + return permutation; - for(int dim=0; dim +PHX::MDField +generateSurfacePermutations(const int num_cells, + const SubcellConnectivity face_connectivity, + PHX::MDField surface_points, + PHX::MDField surface_rotation_matrices) - for(int dim=0; dim<3; ++dim){ - surface_rotation_matrices_k(cell,point,0,dim) = normal[dim]; - surface_rotation_matrices_k(cell,point,1,dim) = transverse[dim]; - surface_rotation_matrices_k(cell,point,2,dim) = binormal[dim]; - } - }); - PHX::Device::execution_space().fence(); - } +{ - // ========================================================= - // Enforce alignment across surface quadrature points + // The challenge for this call is handling wedge-based periodic boundaries + // We need to make sure that we can align points along faces that are rotated with respect to one another. + // Below we will see an algorithm that rotates two boundaries into a shared frame, then figures out + // how to permute the points on one face to line up with the other. - const int num_points = ip_coordinates.extent_int(1); + const int num_points = surface_points.extent_int(1); const int num_faces_per_cell = face_connectivity.numSubcellsOnCellHost(0); const int num_points_per_face = num_points / num_faces_per_cell; + const int cell_dim = surface_points.extent(2); + + MDFieldArrayFactory af("",true); + + // This will store the permutations + auto permutation = af.template buildStaticArray("permutation", num_cells, num_points); + permutation.deep_copy(0); + + // Fill permutations with trivial values (i.e. no permutation - this will get overwritten for some cells) + Kokkos::parallel_for(num_cells,KOKKOS_LAMBDA (const int & cell) { + for(int point=0; point& in_node_ #define PANZER_DOT(a,b) (a[0]*b[0] + a[1]*b[1] + a[2]*b[2]) #define PANZER_CROSS(a,b,c) {c[0] = a[1]*b[2] - a[2]*b[1]; c[1] = a[2]*b[0] - a[0]*b[2]; c[2] = a[0]*b[1] - a[1]*b[0];} - // Reorder index vector. This is scratch space that needs to be - // allocated per thread now. Should use team size instead of num - // faces, but then we need to convert the lambda below to a - // functor to query the policy using the functor. We'll just live - // with the extra memory since it is temporary scratch. - // std::vector point_order(num_points_per_face); - PHX::View point_order("scratch: point_order",face_connectivity.numSubcells(),num_points_per_face); - // Iterate through faces - auto ref_ip_coordinates_k = ref_ip_coordinates.get_view(); - auto ip_coordinates_k = ip_coordinates.get_view(); - auto weighted_measure_k = weighted_measure.get_view(); - auto jac_k = jac.get_view(); - auto jac_det_k = jac_det.get_view(); - auto jac_inv_k = jac_inv.get_view(); - auto surface_normals_k = surface_normals.get_view(); - auto surface_rotation_matrices_k = surface_rotation_matrices.get_view(); Kokkos::parallel_for("face iteration",face_connectivity.numSubcells(),KOKKOS_LAMBDA (const int face) { // Cells for sides 0 and 1 const int cell_0 = face_connectivity.cellForSubcell(face,0); @@ -620,9 +438,9 @@ generateSurfaceCubatureValues(const PHX::MDField& in_node_ for(int face_point=0; face_point& in_node_ const int example_point_1 = lidx_1*num_points_per_face; // Load the transverse and binormal for the 0 cell (default) - Scalar t[3] = {surface_rotation_matrices_k(cell_0,example_point_0,1,0), surface_rotation_matrices_k(cell_0,example_point_0,1,1), surface_rotation_matrices_k(cell_0,example_point_0,1,2)}; - Scalar b[3] = {surface_rotation_matrices_k(cell_0,example_point_0,2,0), surface_rotation_matrices_k(cell_0,example_point_0,2,1), surface_rotation_matrices_k(cell_0,example_point_0,2,2)}; + Scalar t[3] = {surface_rotation_matrices(cell_0,example_point_0,1,0), surface_rotation_matrices(cell_0,example_point_0,1,1), surface_rotation_matrices(cell_0,example_point_0,1,2)}; + Scalar b[3] = {surface_rotation_matrices(cell_0,example_point_0,2,0), surface_rotation_matrices(cell_0,example_point_0,2,1), surface_rotation_matrices(cell_0,example_point_0,2,2)}; // In case the faces are not antiparallel (e.g. periodic wedge), we may need to change the transverse and binormal { // Get the normals for each face for constructing one of the plane vectors (transverse) - const Scalar n0[3] = {surface_rotation_matrices_k(cell_0,example_point_0,0,0), surface_rotation_matrices_k(cell_0,example_point_0,0,1), surface_rotation_matrices_k(cell_0,example_point_0,0,2)}; - const Scalar n1[3] = {surface_rotation_matrices_k(cell_1,example_point_1,0,0), surface_rotation_matrices_k(cell_1,example_point_1,0,1), surface_rotation_matrices_k(cell_1,example_point_1,0,2)}; + const Scalar n0[3] = {surface_rotation_matrices(cell_0,example_point_0,0,0), surface_rotation_matrices(cell_0,example_point_0,0,1), surface_rotation_matrices(cell_0,example_point_0,0,2)}; + const Scalar n1[3] = {surface_rotation_matrices(cell_1,example_point_1,0,0), surface_rotation_matrices(cell_1,example_point_1,0,1), surface_rotation_matrices(cell_1,example_point_1,0,2)}; // n_0*n_1 == -1 (antiparallel), n_0*n_1 == 1 (parallel - bad), |n_0*n_1| < 1 (other) const Scalar n0_dot_n1 = PANZER_DOT(n0,n1); @@ -676,7 +494,7 @@ generateSurfaceCubatureValues(const PHX::MDField& in_node_ t[1] /= mag_t; t[2] /= mag_t; - // The binormal will be the sum of the normals (does not need to be a right handed system) + // The binormal will be the sum of the normals (does not need to be a right handed system) b[0] = n0[0] + n1[0]; b[1] = n0[1] + n1[1]; b[2] = n0[2] + n1[2]; @@ -707,15 +525,12 @@ generateSurfaceCubatureValues(const PHX::MDField& in_node_ // Load position shifted by face center for(int dim=0; dim& in_node_ // Load position shifted by face center for(int dim=0; dim& in_node_ // Find the distance squared between p0 and p1 const Scalar p012 = (p0[0]-p1[0])*(p0[0]-p1[0]) + (p0[1]-p1[1])*(p0[1]-p1[1]); + // TODO: Should this be a constant value, or should we just find the minimum point? // If the distance, compared to the size of the cell, is small, we assume these are the same points if(p012 / r2 < 1.e-12){ - point_order(face,face_point_1) = face_point_0; + permutation(cell_1,lidx_1*num_points_per_face+face_point_0) = point_1; break; } - // Big problem - there wan't a point that aligned properly - KOKKOS_ASSERT(face_point_0 != num_points_per_face-1); - - } - } - - // Now re-order the points on face 1 to correct the alignment - const int point_offset = lidx_1*num_points_per_face; - for( int face_point_1 = 0; face_point_1 < num_points_per_face - 1; ++face_point_1 ){ - // While the point is not yet in place - keep swapping until it is in place (N^2 operations - not great) - while( face_point_1 != point_order(face,face_point_1) ){ - // We need to swap with the component in this position - const int face_point_0 = point_order(face,face_point_1); - panzer::swapQuadraturePoints(cell_1,point_offset+face_point_1,point_offset+face_point_0, - ref_ip_coordinates_k, - ip_coordinates_k, - weighted_measure_k, - jac_k, - jac_det_k, - jac_inv_k, - surface_normals_k, - surface_rotation_matrices_k); - Scalar tmp = point_order(face,face_point_1); - point_order(face,face_point_1) = point_order(face,face_point_0); - point_order(face,face_point_0) = tmp; } } - }); PHX::Device::execution_space().fence(); } @@ -774,472 +564,1240 @@ generateSurfaceCubatureValues(const PHX::MDField& in_node_ #undef PANZER_DOT #undef PANZER_CROSS - // ========================================================= + return permutation; - // I'm not sure if these should exist for surface integrals, but here we go! - - // Shakib contravarient metric tensor - { - auto contravarient_k = this->contravarient.get_static_view(); - auto covarient_k = this->covarient.get_static_view(); - auto jac_k = this->jac.get_static_view(); - Kokkos::MDRangePolicy> policy({0,0},{(int)num_cells,(int)contravarient.extent(1)}); - Kokkos::parallel_for("covarient metric tensor",policy,KOKKOS_LAMBDA (const int cell,const int ip) { - // zero out matrix - for (size_type i = 0; i < contravarient_k.extent(2); ++i) - for (size_type j = 0; j < contravarient_k.extent(3); ++j) - covarient_k(cell,ip,i,j) = 0.0; - - // g^{ij} = \frac{\parital x_i}{\partial \chi_\alpha}\frac{\parital x_j}{\partial \chi_\alpha} - for (size_type i = 0; i < contravarient_k.extent(2); ++i) { - for (size_type j = 0; j < contravarient_k.extent(3); ++j) { - for (size_type alpha = 0; alpha < contravarient_k.extent(2); ++alpha) { - covarient_k(cell,ip,i,j) += jac_k(cell,ip,i,alpha) * jac_k(cell,ip,j,alpha); - } - } - } - }); - PHX::Device::execution_space().fence(); - } +} - { - auto s_contravarient = Kokkos::subview(contravarient.get_view(), std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL,Kokkos::ALL); - auto s_covarient = Kokkos::subview(covarient.get_view(), std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL,Kokkos::ALL); - Intrepid2::RealSpaceTools::inverse(s_contravarient, s_covarient); - PHX::Device::execution_space().fence(); - } +} - // norm of g_ij - { - auto contravarient_k = this->contravarient.get_static_view(); - auto norm_contravarient_k = this->norm_contravarient.get_static_view(); - Kokkos::MDRangePolicy> policy({0,0},{(int)num_cells,(int)contravarient.extent(1)}); - Kokkos::parallel_for("covarient metric tensor",policy,KOKKOS_LAMBDA (const int cell,const int ip) { - norm_contravarient_k(cell,ip) = 0.0; - for (size_type i = 0; i < contravarient_k.extent(2); ++i) { - for (size_type j = 0; j < contravarient_k.extent(3); ++j) { - norm_contravarient_k(cell,ip) += contravarient_k(cell,ip,i,j) * contravarient_k(cell,ip,i,j); - } - } - norm_contravarient_k(cell,ip) = Kokkos::Experimental::sqrt(norm_contravarient_k(cell,ip)); - }); - PHX::Device::execution_space().fence(); - } +//template +//using UnmanagedDynRankView = Kokkos::DynRankView::type,PHX::Device,Kokkos::MemoryTraits>; +template +IntegrationValues2:: +IntegrationValues2(const std::string & pre, + const bool allocArrays): + num_cells_(0), + num_evaluate_cells_(0), + num_virtual_cells_(-1), + requires_permutation_(false), + alloc_arrays_(allocArrays), + prefix_(pre), + ddims_(1,0) +{ + resetArrays(); } +template +void +IntegrationValues2:: +resetArrays() +{ + cub_points_evaluated_ = false; + side_cub_points_evaluated_ = false; + cub_weights_evaluated_ = false; + node_coordinates_evaluated_ = false; + jac_evaluated_ = false; + jac_inv_evaluated_ = false; + jac_det_evaluated_ = false; + weighted_measure_evaluated_ = false; + weighted_normals_evaluated_ = false; + surface_normals_evaluated_ = false; + surface_rotation_matrices_evaluated_ = false; + covarient_evaluated_ = false; + contravarient_evaluated_ = false; + norm_contravarient_evaluated_ = false; + ip_coordinates_evaluated_ = false; + ref_ip_coordinates_evaluated_ = false; + + // TODO: We need to clear the views +} template void IntegrationValues2:: -evaluateRemainingValues(const PHX::MDField& in_node_coordinates, - const int in_num_cells) +setupArrays(const Teuchos::RCP& ir) { - Intrepid2::CellTools cell_tools; + MDFieldArrayFactory af(prefix_,alloc_arrays_); - // copy the dynamic data structures into the static data structures - { - Kokkos::deep_copy(cub_weights.get_static_view(),dyn_cub_weights.get_view()); - Kokkos::deep_copy(cub_points.get_static_view(),dyn_cub_points.get_view()); - } + typedef panzer::IntegrationDescriptor ID; - if (int_rule->isSide()) { - Kokkos::deep_copy(side_cub_points.get_static_view(),dyn_side_cub_points.get_view()); - } + int_rule = ir; - const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells; + const int num_nodes = ir->topology->getNodeCount(); + const int num_cells = ir->workset_size; + const int num_space_dim = ir->topology->getDimension(); - { - size_type num_nodes = in_node_coordinates.extent(1); - size_type num_dims = in_node_coordinates.extent(2); - auto node_coordinates_k = node_coordinates.get_view(); - auto in_node_coordinates_k = in_node_coordinates.get_view(); - - Kokkos::MDRangePolicy> policy({0,0,0},{(int)num_cells,(int)num_nodes,(int)num_dims}); - Kokkos::parallel_for("node coordinates",policy,KOKKOS_LAMBDA (const int cell,const int node,const int dim) { - node_coordinates_k(cell,node,dim) = in_node_coordinates_k(cell,node,dim); - }); - PHX::Device::execution_space().fence(); + // Specialize content if this is quadrature at a node + const bool is_node_rule = (num_space_dim==1) and ir->isSide(); + if(not is_node_rule) { + TEUCHOS_ASSERT(ir->getType() != ID::NONE); + intrepid_cubature = getIntrepidCubature(*ir); } - auto s_in_node_coordinates = Kokkos::subview(in_node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL()); - auto s_jac = Kokkos::subview(jac.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); - cell_tools.setJacobian(jac.get_view(), - cub_points.get_view(), - node_coordinates.get_view(), - *(int_rule->topology)); - PHX::Device::execution_space().fence(); + const int num_ip = is_node_rule ? 1 : ir->num_points; - auto s_jac_inv = Kokkos::subview(jac_inv.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); - cell_tools.setJacobianInv(s_jac_inv, s_jac); + cub_points = af.template buildStaticArray("cub_points",num_ip, num_space_dim); - auto s_jac_det = Kokkos::subview(jac_det.get_view(),std::make_pair(0,num_cells),Kokkos::ALL()); - cell_tools.setJacobianDet(s_jac_det, s_jac); + if (ir->isSide() && ir->cv_type == "none") + side_cub_points = af.template buildStaticArray("side_cub_points",num_ip,ir->side_topology->getDimension()); - PHX::Device::execution_space().fence(); + cub_weights = af.template buildStaticArray("cub_weights",num_ip); - auto s_weighted_measure = Kokkos::subview(weighted_measure.get_view(),std::make_pair(0,num_cells),Kokkos::ALL()); - if (!int_rule->isSide()) { - Intrepid2::FunctionSpaceTools:: - computeCellMeasure(s_weighted_measure, s_jac_det, cub_weights.get_view()); - } - else if(int_rule->spatial_dimension==3) { - Intrepid2::FunctionSpaceTools:: - computeFaceMeasure(s_weighted_measure, s_jac, cub_weights.get_view(), - int_rule->side, *int_rule->topology, - scratch_for_compute_side_measure.get_view()); - } - else if(int_rule->spatial_dimension==2) { - Intrepid2::FunctionSpaceTools:: - computeEdgeMeasure(s_weighted_measure, s_jac, cub_weights.get_view(), - int_rule->side,*int_rule->topology, - scratch_for_compute_side_measure.get_view()); - } - else TEUCHOS_ASSERT(false); + node_coordinates = af.template buildStaticArray("node_coordinates",num_cells, num_nodes, num_space_dim); - PHX::Device::execution_space().fence(); + jac = af.template buildStaticArray("jac",num_cells, num_ip, num_space_dim,num_space_dim); - // Shakib contravarient metric tensor - { - auto covarient_k = covarient.get_view(); - auto contravarient_k = contravarient.get_view(); - auto jac_k = jac.get_view(); - Kokkos::MDRangePolicy> policy({0,0},{num_cells,static_cast(contravarient.extent(1))}); - Kokkos::parallel_for("evalaute covarient metric tensor",policy,KOKKOS_LAMBDA (const int cell,const int ip) { - // zero out matrix - for (int i = 0; i < static_cast(contravarient_k.extent(2)); ++i) - for (int j = 0; j < static_cast(contravarient_k.extent(3)); ++j) - covarient_k(cell,ip,i,j) = 0.0; - - // g^{ij} = \frac{\parital x_i}{\partial \chi_\alpha}\frac{\parital x_j}{\partial \chi_\alpha} - for (int i = 0; i < static_cast(contravarient_k.extent(2)); ++i) { - for (int j = 0; j < static_cast(contravarient_k.extent(3)); ++j) { - for (int alpha = 0; alpha < static_cast(contravarient_k.extent(2)); ++alpha) { - covarient_k(cell,ip,i,j) += jac_k(cell,ip,i,alpha) * jac_k(cell,ip,j,alpha); - } - } - } - }); - PHX::Device::execution_space().fence(); - } + jac_inv = af.template buildStaticArray("jac_inv",num_cells, num_ip, num_space_dim,num_space_dim); - auto s_covarient = Kokkos::subview(covarient.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); - auto s_contravarient = Kokkos::subview(contravarient.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); - Intrepid2::RealSpaceTools::inverse(s_contravarient, s_covarient); - PHX::Device::execution_space().fence(); + jac_det = af.template buildStaticArray("jac_det",num_cells, num_ip); + + weighted_measure = af.template buildStaticArray("weighted_measure",num_cells, num_ip); + + covarient = af.template buildStaticArray("covarient",num_cells, num_ip, num_space_dim,num_space_dim); + + contravarient = af.template buildStaticArray("contravarient",num_cells, num_ip, num_space_dim,num_space_dim); + + norm_contravarient = af.template buildStaticArray("norm_contravarient",num_cells, num_ip); + + ip_coordinates = af.template buildStaticArray("ip_coordiantes",num_cells, num_ip,num_space_dim); + + ref_ip_coordinates = af.template buildStaticArray("ref_ip_coordinates",num_cells, num_ip,num_space_dim); + + weighted_normals = af.template buildStaticArray("weighted_normal",num_cells,num_ip,num_space_dim); + + surface_normals = af.template buildStaticArray("surface_normals",num_cells, num_ip,num_space_dim); + + surface_rotation_matrices = af.template buildStaticArray("surface_rotation_matrices",num_cells, num_ip,3,3); - // norm of g_ij - { - auto contravarient_k = contravarient.get_view(); - auto norm_contravarient_k = norm_contravarient.get_view(); - Kokkos::MDRangePolicy> policy({0,0},{num_cells,static_cast(contravarient.extent(1))}); - Kokkos::parallel_for("evaluate norm_contravarient",policy,KOKKOS_LAMBDA (const int cell,const int ip) { - norm_contravarient_k(cell,ip) = 0.0; - for (int i = 0; i < static_cast(contravarient_k.extent(2)); ++i) { - for (int j = 0; j < static_cast(contravarient_k.extent(3)); ++j) { - norm_contravarient_k(cell,ip) += contravarient_k(cell,ip,i,j) * contravarient_k(cell,ip,i,j); - } - } - norm_contravarient_k(cell,ip) = Kokkos::Experimental::sqrt(norm_contravarient_k(cell,ip)); - }); - PHX::Device::execution_space().fence(); - } } -// Find the permutation that maps the set of points coords to other_coords. To -// avoid possible finite precision issues, == is not used, but rather -// min(norm(.)). + +// *********************************************************** +// * Evaluation of values - NOT specialized +// *********************************************************** template -static void -permuteToOther(const PHX::MDField& coords, - const PHX::MDField& other_coords, - std::vector >::size_type>& permutation) +void IntegrationValues2:: +evaluateValues(const PHX::MDField & in_node_coordinates, + const int in_num_cells, + const Teuchos::RCP & face_connectivity, + const int num_virtual_cells) { - typedef typename ArrayTraits >::size_type size_type; - // We can safely assume: (1) The permutation is the same for every cell in - // the workset. (2) The first workset has valid data. Hence we operate only - // on cell 0. - const size_type cell = 0; - const size_type num_ip = coords.extent(1), num_dim = coords.extent(2); - permutation.resize(num_ip); - std::vector taken(num_ip, 0); - - auto coords_view = coords.get_view(); - auto coords_h = Kokkos::create_mirror_view(coords_view); - Kokkos::deep_copy(coords_h, coords_view); - - auto other_coords_view = other_coords.get_view(); - auto other_coords_h = Kokkos::create_mirror_view(other_coords_view); - Kokkos::deep_copy(other_coords_h, other_coords_view); - - for (size_type ip = 0; ip < num_ip; ++ip) { - // Find an other point to associate with ip. - size_type i_min = 0; - Scalar d_min = -1; - for (size_type other_ip = 0; other_ip < num_ip; ++other_ip) { - // For speed, skip other points that are already associated. - if (taken[other_ip]) continue; - // Compute the distance between the two points. - Scalar d(0); - for (size_type dim = 0; dim < num_dim; ++dim) { - const Scalar diff = coords_h(cell, ip, dim) - other_coords_h(cell, other_ip, dim); - d += diff*diff; - } - if (d_min < 0 || d < d_min) { - d_min = d; - i_min = other_ip; - } - } - // Record the match. - permutation[ip] = i_min; - // This point is now taken. - taken[i_min] = 1; - } + + setup(int_rule, in_node_coordinates, in_num_cells); + + // Setup permutations (only required for surface integrators) + using ID=panzer::IntegrationDescriptor; + if((int_rule->getType() == ID::SURFACE) and (face_connectivity != Teuchos::null)) + setupPermutations(face_connectivity, num_virtual_cells); + + // Evaluate everything once permutations are generated + evaluateEverything(); + } template -void IntegrationValues2:: +void +IntegrationValues2:: evaluateValues(const PHX::MDField& in_node_coordinates, const PHX::MDField& other_ip_coordinates, const int in_num_cells) { - const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells; - if (int_rule->cv_type == "none") { + setup(int_rule, in_node_coordinates, in_num_cells); - getCubature(in_node_coordinates, in_num_cells); + // Setup permutations + setupPermutations(other_ip_coordinates); - { - // Determine the permutation. - std::vector permutation(other_ip_coordinates.extent(1)); - permuteToOther(ip_coordinates, other_ip_coordinates, permutation); - // Apply the permutation to the cubature arrays. - MDFieldArrayFactory af(prefix, alloc_arrays); - { - const size_type num_ip = dyn_cub_points.extent(0); - { - const size_type num_dim = dyn_side_cub_points.extent(1); - DblArrayDynamic old_dyn_side_cub_points = af.template buildArray( - "old_dyn_side_cub_points", num_ip, num_dim); - old_dyn_side_cub_points.deep_copy(dyn_side_cub_points); - - auto dyn_side_cub_points_h = Kokkos::create_mirror_view(as_view(dyn_side_cub_points)); - auto old_dyn_side_cub_points_h = Kokkos::create_mirror_view(as_view(old_dyn_side_cub_points)); - Kokkos::deep_copy(dyn_side_cub_points_h, as_view(dyn_side_cub_points)); - Kokkos::deep_copy(old_dyn_side_cub_points_h, as_view(old_dyn_side_cub_points)); - - for (size_type ip = 0; ip < num_ip; ++ip) - if (ip != permutation[ip]) - for (size_type dim = 0; dim < num_dim; ++dim) - dyn_side_cub_points_h(ip, dim) = old_dyn_side_cub_points_h(permutation[ip], dim); - - Kokkos::deep_copy(as_view(dyn_side_cub_points), dyn_side_cub_points_h); - } - { - const size_type num_dim = dyn_cub_points.extent(1); - DblArrayDynamic old_dyn_cub_points = af.template buildArray( - "old_dyn_cub_points", num_ip, num_dim); - old_dyn_cub_points.deep_copy(dyn_cub_points); - - auto dyn_cub_points_h = Kokkos::create_mirror_view(as_view(dyn_cub_points)); - auto old_dyn_cub_points_h = Kokkos::create_mirror_view(as_view(old_dyn_cub_points)); - Kokkos::deep_copy(dyn_cub_points_h, as_view(dyn_cub_points)); - Kokkos::deep_copy(old_dyn_cub_points_h, as_view(old_dyn_cub_points)); - - for (size_type ip = 0; ip < num_ip; ++ip) - if (ip != permutation[ip]) - for (size_type dim = 0; dim < num_dim; ++dim) - dyn_cub_points_h(ip, dim) = old_dyn_cub_points_h(permutation[ip], dim); - - Kokkos::deep_copy(as_view(dyn_cub_points), dyn_cub_points_h); - } - { - DblArrayDynamic old_dyn_cub_weights = af.template buildArray( - "old_dyn_cub_weights", num_ip); - old_dyn_cub_weights.deep_copy(dyn_cub_weights); + // Evaluate everything once permutations are generated + evaluateEverything(); - auto dyn_cub_weights_h = Kokkos::create_mirror_view(as_view(dyn_cub_weights)); - auto old_dyn_cub_weights_h = Kokkos::create_mirror_view(as_view(old_dyn_cub_weights)); - Kokkos::deep_copy(dyn_cub_weights_h, as_view(dyn_cub_weights)); - Kokkos::deep_copy(old_dyn_cub_weights_h, as_view(old_dyn_cub_weights)); +} - for (size_type ip = 0; ip < dyn_cub_weights.extent(0); ++ip) - if (ip != permutation[ip]) - dyn_cub_weights_h(ip) = old_dyn_cub_weights_h(permutation[ip]); +template +void +IntegrationValues2:: +setupPermutations(const Teuchos::RCP & face_connectivity, + const int num_virtual_cells) +{ + TEUCHOS_ASSERT(not int_rule->isSide()); + TEUCHOS_ASSERT(face_connectivity != Teuchos::null); + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() != panzer::IntegrationDescriptor::SURFACE, + "IntegrationValues2::setupPermutations : Face connectivity based permutations are only required for surface integration schemes"); + TEUCHOS_TEST_FOR_EXCEPT_MSG(num_virtual_cells_ >= 0, + "IntegrationValues2::setupPermutations : Number of virtual cells for surface permuations must be >= 0"); + resetArrays(); + side_connectivity_ = face_connectivity; + num_virtual_cells_ = num_virtual_cells; + requires_permutation_ = false; + permutations_ = generateSurfacePermutations(num_evaluate_cells_,*face_connectivity, getCubaturePoints(false,true), getSurfaceRotationMatrices(false,true)); + requires_permutation_ = true; +} - Kokkos::deep_copy(as_view(dyn_cub_weights), old_dyn_cub_weights_h); - } - } - { - const size_type num_ip = ip_coordinates.extent(1); - const size_type num_dim = ip_coordinates.extent(2); - Array_CellIPDim old_ip_coordinates = af.template buildStaticArray( - "old_ip_coordinates", ip_coordinates.extent(0), num_ip, num_dim); - - auto ip_coordinates_h = Kokkos::create_mirror_view(as_view(ip_coordinates)); - auto old_ip_coordinates_h = Kokkos::create_mirror_view(as_view(old_ip_coordinates)); - Kokkos::deep_copy(old_ip_coordinates_h, as_view(ip_coordinates)); - Kokkos::deep_copy(ip_coordinates_h, as_view(ip_coordinates)); - - for (int cell = 0; cell < num_cells; ++cell) - for (size_type ip = 0; ip < num_ip; ++ip) - if (ip != permutation[ip]) - for (size_type dim = 0; dim < num_dim; ++dim) - ip_coordinates_h(cell, ip, dim) = old_ip_coordinates_h(cell, permutation[ip], dim); - - Kokkos::deep_copy(as_view(ip_coordinates), ip_coordinates_h); - } - // All subsequent calculations inherit the permutation. - } +template +void +IntegrationValues2:: +setupPermutations(const PHX::MDField & other_ip_coordinates) +{ + resetArrays(); + requires_permutation_ = false; + permutations_ = generatePermutations(num_evaluate_cells_, getCubaturePoints(false,true), other_ip_coordinates); + requires_permutation_ = true; +} + + +template +void +IntegrationValues2:: +setup(const Teuchos::RCP& ir, + const PHX::MDField & vertex_coordinates, + const int num_cells) +{ + + // Clear arrays just in case we are rebuilding this object + resetArrays(); + + num_cells_ = vertex_coordinates.extent(0); + num_evaluate_cells_ = num_cells < 0 ? vertex_coordinates.extent(0) : num_cells; + int_rule = ir; + + TEUCHOS_ASSERT(ir->getType() != IntegrationDescriptor::NONE); + intrepid_cubature = getIntrepidCubature(*ir); + + // Copy node coordinates into owned allocation + { + MDFieldArrayFactory af(prefix_,true); + + const int num_space_dim = int_rule->topology->getDimension(); + const int num_vertexes = int_rule->topology->getNodeCount(); + TEUCHOS_ASSERT(static_cast(vertex_coordinates.extent(1)) == num_vertexes); + + auto aux = af.template buildStaticArray("node_coordinates",num_cells_,num_vertexes, num_space_dim); + Kokkos::MDRangePolicy> policy({0,0,0},{num_evaluate_cells_,num_vertexes,num_space_dim}); + Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int & cell, const int & point, const int & dim){ + aux(cell,point,dim) = vertex_coordinates(cell,point,dim); + }); + PHX::Device::execution_space().fence(); + node_coordinates = aux; + } + +} + +template +typename IntegrationValues2::ConstArray_IPDim +IntegrationValues2:: +getUniformCubaturePointsRef(const bool cache, + const bool force, + const bool apply_permutation) const +{ - evaluateRemainingValues(in_node_coordinates, in_num_cells); + if(cub_points_evaluated_ and not force) + return cub_points; + + Intrepid2::CellTools cell_tools; + MDFieldArrayFactory af(prefix_,true); + + int num_space_dim = int_rule->topology->getDimension(); + int num_ip = int_rule->num_points; + + auto aux = af.template buildStaticArray("cub_points",num_ip, num_space_dim); + + if (int_rule->isSide() && num_space_dim==1) { + std::cout << "WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do " + << "non-natural integration rules."; + return aux; } + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type != "none", + "IntegrationValues2::getUniformCubaturePointsRef : Cannot build reference cubature points for control volume integration scheme."); + + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() == IntegrationDescriptor::SURFACE, + "IntegrationValues2::getUniformCubaturePointsRef : Cannot build reference cubature points for surface integration scheme."); + + auto weights = af.template buildStaticArray("cub_weights",num_ip); + + if (!int_rule->isSide()) + intrepid_cubature->getCubature(aux.get_view(), weights.get_view()); else { + auto s_cub_points = af.template buildStaticArray("side_cub_points",num_ip, num_space_dim-1); - getCubatureCV(in_node_coordinates, in_num_cells); - - // Determine the permutation. - std::vector permutation(other_ip_coordinates.extent(1)); - permuteToOther(ip_coordinates, other_ip_coordinates, permutation); - - // Apply the permutation to the cubature arrays. - MDFieldArrayFactory af(prefix, alloc_arrays); - { - const size_type workset_size = ip_coordinates.extent(0), num_ip = ip_coordinates.extent(1), - num_dim = ip_coordinates.extent(2); - Array_CellIPDim old_ip_coordinates = af.template buildStaticArray( - "old_ip_coordinates", workset_size, num_ip, num_dim); - Kokkos::deep_copy(old_ip_coordinates.get_static_view(), ip_coordinates.get_static_view()); - Array_CellIPDim old_weighted_normals = af.template buildStaticArray( - "old_weighted_normals", workset_size, num_ip, num_dim); - Array_CellIP old_weighted_measure = af.template buildStaticArray( - "old_weighted_measure", workset_size, num_ip); - if (int_rule->cv_type == "side") - Kokkos::deep_copy(old_weighted_normals.get_static_view(), weighted_normals.get_static_view()); - else - Kokkos::deep_copy(old_weighted_measure.get_static_view(), weighted_measure.get_static_view()); - for (int cell = 0; cell < num_cells; ++cell) - { - for (size_type ip = 0; ip < num_ip; ++ip) - { - if (ip != permutation[ip]) { - if (int_rule->cv_type == "boundary" || int_rule->cv_type == "volume") - weighted_measure(cell, ip) = old_weighted_measure(cell, permutation[ip]); - for (size_type dim = 0; dim < num_dim; ++dim) - { - ip_coordinates(cell, ip, dim) = old_ip_coordinates(cell, permutation[ip], dim); - if (int_rule->cv_type == "side") - weighted_normals(cell, ip, dim) = old_weighted_normals(cell, permutation[ip], dim); - - } - } - } - } - } + intrepid_cubature->getCubature(s_cub_points.get_view(), weights.get_view()); + cell_tools.mapToReferenceSubcell(aux.get_view(), s_cub_points.get_view(), num_space_dim-1, int_rule->getSide(), *(int_rule->topology)); + } + + PHX::Device::execution_space().fence(); + + if(apply_permutation and requires_permutation_) + applyBasePermutation(aux, permutations_); - evaluateValuesCV(in_node_coordinates, in_num_cells); + if(cache){ + cub_points = aux; + cub_points_evaluated_ = true; } + + return aux; + } template -void IntegrationValues2:: -getCubatureCV(const PHX::MDField& in_node_coordinates, - const int in_num_cells) +typename IntegrationValues2::ConstArray_IPDim +IntegrationValues2:: +getUniformSideCubaturePointsRef(const bool cache, + const bool force, + const bool apply_permutation) const { + + if(side_cub_points_evaluated_ and not force) + return side_cub_points; + + Intrepid2::CellTools cell_tools; + MDFieldArrayFactory af(prefix_,true); + int num_space_dim = int_rule->topology->getDimension(); + int num_ip = int_rule->num_points; + + auto aux = af.template buildStaticArray("side_cub_points",num_ip, num_space_dim-1); + if (int_rule->isSide() && num_space_dim==1) { std::cout << "WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do " << "non-natural integration rules."; - return; + return aux; } - size_type num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : (size_type) in_num_cells; - std::pair cell_range(0,num_cells); - { - size_type num_nodes = in_node_coordinates.extent(1); - size_type num_dims = in_node_coordinates.extent(2); - auto node_coordinates_k = node_coordinates.get_view(); - auto dyn_node_coordinates_k = dyn_node_coordinates.get_view(); - auto in_node_coordinates_k = in_node_coordinates.get_view(); - - Kokkos::MDRangePolicy> policy({0,0,0},{num_cells,num_nodes,num_dims}); - Kokkos::parallel_for("getCubatureCV: node coordinates",policy,KOKKOS_LAMBDA (const int cell,const int node,const int dim) { - node_coordinates_k(cell,node,dim) = in_node_coordinates_k(cell,node,dim); - dyn_node_coordinates_k(cell,node,dim) = Sacado::ScalarValue::eval(in_node_coordinates_k(cell,node,dim)); - }); - PHX::Device::execution_space().fence(); - } + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type != "none", + "IntegrationValues2::getUniformSideCubaturePointsRef : Cannot build reference cubature points for control volume integration scheme."); + + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() == IntegrationDescriptor::SURFACE, + "IntegrationValues2::getUniformSideCubaturePointsRef : Cannot build reference cubature points for surface integration scheme."); + + TEUCHOS_TEST_FOR_EXCEPT_MSG(!int_rule->isSide(), + "IntegrationValues2::getUniformSideCubaturePointsRef : Requested side points, which is not supported by integration rule."); + + auto weights = af.template buildStaticArray("cub_weights",num_ip); + + intrepid_cubature->getCubature(aux.get_view(), weights.get_view()); - auto s_dyn_phys_cub_points = Kokkos::subdynrankview(dyn_phys_cub_points.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); - auto s_dyn_node_coordinates = Kokkos::subdynrankview(dyn_node_coordinates.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); - if (int_rule->cv_type == "side") { - auto s_dyn_phys_cub_norms = Kokkos::subdynrankview(dyn_phys_cub_norms.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); - intrepid_cubature->getCubature(s_dyn_phys_cub_points,s_dyn_phys_cub_norms,s_dyn_node_coordinates); + PHX::Device::execution_space().fence(); + + if(apply_permutation and requires_permutation_) + applyBasePermutation(aux, permutations_); + + if(cache){ + side_cub_points = aux; + side_cub_points_evaluated_ = true; } - else { - auto s_dyn_phys_cub_weights = Kokkos::subdynrankview(dyn_phys_cub_weights.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); - intrepid_cubature->getCubature(s_dyn_phys_cub_points,s_dyn_phys_cub_weights,s_dyn_node_coordinates); + + return aux; +} + +template +typename IntegrationValues2::ConstArray_IP +IntegrationValues2:: +getUniformCubatureWeightsRef(const bool cache, + const bool force, + const bool apply_permutation) const +{ + + if(cub_weights_evaluated_ and not force) + return cub_weights; + + Intrepid2::CellTools cell_tools; + MDFieldArrayFactory af(prefix_,true); + + int num_space_dim = int_rule->topology->getDimension(); + int num_ip = int_rule->num_points; + + auto aux = af.template buildStaticArray("cub_weights",num_ip); + + if (int_rule->isSide() && num_space_dim==1) { + std::cout << "WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do " + << "non-natural integration rules."; + return aux; } - // size_type num_cells = dyn_phys_cub_points.extent(0); - size_type num_ip =dyn_phys_cub_points.extent(1); - size_type num_dims = dyn_phys_cub_points.extent(2); - auto weighted_measure_k = weighted_measure.get_view(); - auto dyn_phys_cub_weights_k = dyn_phys_cub_weights.get_view(); - auto ip_coordinates_k = ip_coordinates.get_view(); - auto dyn_phys_cub_points_k = dyn_phys_cub_points.get_view(); - auto weighted_normals_k = weighted_normals.get_view(); - auto dyn_phys_cub_norms_k = dyn_phys_cub_norms.get_view(); - bool is_side = false; - if (int_rule->cv_type == "side") - is_side = true; - - Kokkos::MDRangePolicy> policy({0,0},{num_cells,num_ip}); - Kokkos::parallel_for("getCubatureCV: weighted measure",policy,KOKKOS_LAMBDA (const int cell,const int ip) { - if (!is_side) - weighted_measure_k(cell,ip) = dyn_phys_cub_weights_k(cell,ip); - for (size_type dim = 0; dim < num_dims; ++dim) { - ip_coordinates_k(cell,ip,dim) = dyn_phys_cub_points_k(cell,ip,dim); - if (is_side) - weighted_normals_k(cell,ip,dim) = dyn_phys_cub_norms_k(cell,ip,dim); - } - }); + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type != "none", + "IntegrationValues2::getUniformCubatureWeightsRef : Cannot build reference cubature weights for control volume integration scheme."); + + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() == IntegrationDescriptor::SURFACE, + "IntegrationValues2::getUniformCubatureWeightsRef : Cannot build reference cubature weights for surface integration scheme."); + + auto points = af.template buildStaticArray("cub_points",num_ip,num_space_dim); + + intrepid_cubature->getCubature(points.get_view(), aux.get_view()); + PHX::Device::execution_space().fence(); + + if(apply_permutation and requires_permutation_) + applyBasePermutation(aux, permutations_); + + if(cache){ + cub_weights = aux; + cub_weights_evaluated_ = true; + } + + return aux; + } template -void IntegrationValues2:: -evaluateValuesCV(const PHX::MDField& in_node_coordinates, - const int in_num_cells) +typename IntegrationValues2::ConstArray_CellBASISDim +IntegrationValues2:: +getNodeCoordinates() const +{ + return node_coordinates; +} + +template +typename IntegrationValues2::ConstArray_CellIPDimDim +IntegrationValues2:: +getJacobian(const bool cache, + const bool force) const { + if(jac_evaluated_ and not force) + return jac; Intrepid2::CellTools cell_tools; + MDFieldArrayFactory af(prefix_,true); + + int num_space_dim = int_rule->topology->getDimension(); + int num_ip = int_rule->num_points; + + // Don't forget that since we are not caching this, we have to make sure the managed view remains alive while we use the non-const wrapper + auto const_ref_coord = getCubaturePointsRef(false,force); + auto ref_coord = PHX::getNonConstDynRankViewFromConstMDField(const_ref_coord); + auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates()); + auto aux = af.template buildStaticArray("jac",num_cells_, num_ip, num_space_dim,num_space_dim); + + const auto cell_range = std::make_pair(0,num_evaluate_cells_); + auto s_ref_coord = Kokkos::subview(ref_coord, cell_range,Kokkos::ALL(),Kokkos::ALL()); + auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL()); + auto s_jac = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); + + cell_tools.setJacobian(s_jac, s_ref_coord, s_node_coord,*(int_rule->topology)); + + PHX::Device::execution_space().fence(); + + if(cache){ + jac = aux; + jac_evaluated_ = true; + } + + return aux; - size_type num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : (size_type) in_num_cells; +} + +template +typename IntegrationValues2::ConstArray_CellIPDimDim +IntegrationValues2:: +getJacobianInverse(const bool cache, + const bool force) const +{ - auto s_ref_ip_coordinates = Kokkos::subview(ref_ip_coordinates.get_view(),std::make_pair(0,(int)num_cells),Kokkos::ALL(),Kokkos::ALL()); - auto s_ip_coordinates = Kokkos::subview(ip_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL()); - auto s_node_coordinates = Kokkos::subview(node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL()); + if(jac_inv_evaluated_ and not force) + return jac_inv; - cell_tools.mapToReferenceFrame(s_ref_ip_coordinates, - s_ip_coordinates, - s_node_coordinates, - *(int_rule->topology)); + Intrepid2::CellTools cell_tools; + MDFieldArrayFactory af(prefix_,true); - auto s_jac = Kokkos::subview(jac.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); + const int num_space_dim = int_rule->topology->getDimension(); + const int num_ip = int_rule->num_points; - cell_tools.setJacobian(s_jac, - s_ref_ip_coordinates, - s_node_coordinates, - *(int_rule->topology)); + auto jacobian = getJacobian(false,force); + auto aux = af.template buildStaticArray("jac_inv",num_cells_, num_ip, num_space_dim,num_space_dim); - auto s_jac_inv = Kokkos::subview(jac_inv.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); + const auto cell_range = std::make_pair(0,num_evaluate_cells_); + auto s_jac = Kokkos::subview(jacobian.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); + auto s_jac_inv = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); cell_tools.setJacobianInv(s_jac_inv, s_jac); - auto s_jac_det = Kokkos::subview(jac_det.get_view(),std::make_pair(0,num_cells),Kokkos::ALL()); + PHX::Device::execution_space().fence(); + + if(cache){ + jac_inv = aux; + jac_inv_evaluated_ = true; + } + + return aux; + +} + +template +typename IntegrationValues2::ConstArray_CellIP +IntegrationValues2:: +getJacobianDeterminant(const bool cache, + const bool force) const +{ + + if(jac_det_evaluated_ and not force) + return jac_det; + + Intrepid2::CellTools cell_tools; + MDFieldArrayFactory af(prefix_,true); + + const int num_ip = int_rule->num_points; + + auto jacobian = getJacobian(false,force); + auto aux = af.template buildStaticArray("jac_det",num_cells_, num_ip); + + const auto cell_range = std::make_pair(0,num_evaluate_cells_); + auto s_jac = Kokkos::subview(jacobian.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); + auto s_jac_det = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL()); + + cell_tools.setJacobianDet(s_jac_det, s_jac); + + PHX::Device::execution_space().fence(); + + if(cache){ + jac_det = aux; + jac_det_evaluated_ = true; + } + + return aux; + +} + +template +typename IntegrationValues2::ConstArray_CellIP +IntegrationValues2:: +getWeightedMeasure(const bool cache, + const bool force) const +{ + + if(weighted_measure_evaluated_ and not force) + return weighted_measure; + + Intrepid2::CellTools cell_tools; + MDFieldArrayFactory af(prefix_,true); + + const int num_space_dim = int_rule->topology->getDimension(); + const int num_ip = int_rule->num_points; + + auto aux = af.template buildStaticArray("weighted_measure",num_cells_, num_ip); + + if(int_rule->cv_type != "none"){ + + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type == "side", + "IntegrationValues2::getWeightedMeasure : Weighted measure not available for side control volume integrators. Use getWeightedNormals instead."); + + // CV integration uses a single call to map from physical space to the weighted measure - I assume this is slower than what we do with non-cv integration methods + + auto s_cub_points = af.template buildStaticArray("cub_points",num_evaluate_cells_,num_ip,num_space_dim); + + auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates()); + + const auto cell_range = std::make_pair(0,num_evaluate_cells_); + auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL()); + auto s_weighted_measure = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL()); + + intrepid_cubature->getCubature(s_cub_points.get_view(),s_weighted_measure,s_node_coord); + + } else if(int_rule->getType() == IntegrationDescriptor::SURFACE){ + + const auto & cell_topology = *int_rule->topology; + const int cell_dim = cell_topology.getDimension(); + const int num_sides = (cell_dim==1) ? 2 : cell_topology.getSideCount(); + const int cubature_order = int_rule->order(); + const int num_points_on_side = num_ip / num_sides; + + Intrepid2::DefaultCubatureFactory cubature_factory; + auto jacobian = getJacobian(false,force); + + for(int side=0; side side_cub_weights; + if(cell_dim==1){ + side_cub_weights = Kokkos::DynRankView("side_cub_weights",num_points_on_side); + auto tmp_side_cub_weights_host = Kokkos::create_mirror_view(side_cub_weights); + tmp_side_cub_weights_host(0)=1.; + Kokkos::deep_copy(side_cub_weights,tmp_side_cub_weights_host); + } else { + + // Get the face topology from the cell topology + const shards::CellTopology face_topology(cell_topology.getCellTopologyData(cell_dim-1,side)); + + auto ic = cubature_factory.create(face_topology,cubature_order); + + side_cub_weights = Kokkos::DynRankView("side_cub_weights",num_points_on_side); + auto subcell_cub_points = Kokkos::DynRankView("subcell_cub_points",num_points_on_side,cell_dim-1); + + // Get the reference face points + ic->getCubature(subcell_cub_points, side_cub_weights); + } + + PHX::Device::execution_space().fence(); + + // Iterating over face points + Kokkos::MDRangePolicy> policy({0,0},{num_evaluate_cells_,num_points_on_side}); + + // Calculate measures (quadrature weights in physical space) for this side + auto side_weighted_measure = Kokkos::DynRankView("side_weighted_measure",num_evaluate_cells_,num_points_on_side); + if(cell_dim == 1){ + Kokkos::deep_copy(side_weighted_measure, side_cub_weights(0)); + } else { + + // Copy from complete jacobian to side jacobian + auto side_jacobian = Kokkos::DynRankView("side_jac",num_evaluate_cells_,num_points_on_side,cell_dim,cell_dim); + + Kokkos::parallel_for("Copy jacobian to side jacobian",policy,KOKKOS_LAMBDA (const int cell,const int point) { + for(int dim=0;dim("scratch_for_compute_measure", num_evaluate_cells_*num_points_on_side*num_space_dim*num_space_dim); + + if(cell_dim == 2){ + Intrepid2::FunctionSpaceTools:: + computeEdgeMeasure(side_weighted_measure, side_jacobian, side_cub_weights, + side,cell_topology, + scratch.get_view()); + PHX::Device::execution_space().fence(); + } else if(cell_dim == 3){ + Intrepid2::FunctionSpaceTools:: + computeFaceMeasure(side_weighted_measure, side_jacobian, side_cub_weights, + side,cell_topology, + scratch.get_view()); + PHX::Device::execution_space().fence(); + } + } + + + // Copy to the main array + Kokkos::parallel_for("copy surface weighted measure values",policy,KOKKOS_LAMBDA (const int cell,const int point) { + aux(cell,point_offset + point) = side_weighted_measure(cell,point); + }); + PHX::Device::execution_space().fence(); + } + + } else { + + auto cell_range = std::make_pair(0,num_evaluate_cells_); + auto s_weighted_measure = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL()); + auto cubature_weights = getUniformCubatureWeightsRef(false,force,false); + + if (!int_rule->isSide()) { + + auto s_jac_det = Kokkos::subview(getJacobianDeterminant(false,force).get_view(),cell_range,Kokkos::ALL()); + Intrepid2::FunctionSpaceTools:: + computeCellMeasure(s_weighted_measure, s_jac_det, cubature_weights.get_view()); + + } else if(int_rule->spatial_dimension==3) { + + auto s_jac = Kokkos::subview(getJacobian(false,force).get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); + auto scratch = af.template buildStaticArray("scratch_for_compute_measure", num_evaluate_cells_*num_ip*num_space_dim*num_space_dim); + Intrepid2::FunctionSpaceTools:: + computeFaceMeasure(s_weighted_measure, s_jac, cubature_weights.get_view(), + int_rule->side, *int_rule->topology, + scratch.get_view()); + + } else if(int_rule->spatial_dimension==2) { + + auto s_jac = Kokkos::subview(getJacobian(false,force).get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()); + auto scratch = af.template buildStaticArray("scratch_for_compute_measure", num_evaluate_cells_*num_ip*num_space_dim*num_space_dim); + Intrepid2::FunctionSpaceTools:: + computeEdgeMeasure(s_weighted_measure, s_jac, cubature_weights.get_view(), + int_rule->side,*int_rule->topology, + scratch.get_view()); + + } else { + TEUCHOS_ASSERT(false); + } + + } + + PHX::Device::execution_space().fence(); + + // Apply permutations if necessary + if(requires_permutation_) + applyPermutation(aux, permutations_); + + if(cache){ + weighted_measure = aux; + weighted_measure_evaluated_ = true; + } + + return aux; + +} + +template +typename IntegrationValues2::ConstArray_CellIPDim +IntegrationValues2:: +getWeightedNormals(const bool cache, + const bool force) const +{ + + if(weighted_normals_evaluated_ and not force) + return weighted_normals; + + Intrepid2::CellTools cell_tools; + MDFieldArrayFactory af(prefix_,true); + + const int num_space_dim = int_rule->topology->getDimension(); + const int num_ip = int_rule->num_points; + + auto aux = af.template buildStaticArray("weighted_normals",num_cells_,num_ip,num_space_dim); + + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() == IntegrationDescriptor::SURFACE, + "IntegrationValues2::getWeightedNormals : Cannot build reference weighted normals for surface integration scheme."); + + auto points = af.template buildStaticArray("cub_points",num_cells_,num_ip,num_space_dim); + + auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates()); + + const auto cell_range = std::make_pair(0,num_evaluate_cells_); + auto s_cub_points = Kokkos::subview(points.get_view(),cell_range, Kokkos::ALL(), Kokkos::ALL()); + auto s_weighted_normals = Kokkos::subview(aux.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL()); + auto s_node_coord = Kokkos::subview(node_coord, cell_range, Kokkos::ALL(), Kokkos::ALL()); + + intrepid_cubature->getCubature(s_cub_points,s_weighted_normals,s_node_coord); + + PHX::Device::execution_space().fence(); + + // Apply permutations if necessary + if(requires_permutation_) + applyPermutation(aux, permutations_); + + if(cache){ + weighted_normals = aux; + weighted_normals_evaluated_ = true; + } + + return aux; + +} + +template +typename IntegrationValues2::ConstArray_CellIPDim +IntegrationValues2:: +getSurfaceNormals(const bool cache, + const bool force) const +{ + + if(surface_normals_evaluated_ and not force) + return surface_normals; + + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->isSide(), + "IntegrationValues2::getSurfaceNormals : This call doesn't work with sides (only surfaces)."); + + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type != "none", + "IntegrationValues2::getSurfaceNormals : This call does not support control volume integration schemes."); + + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() != IntegrationDescriptor::SURFACE, + "IntegrationValues2::getSurfaceNormals : Can only build for surface integrators."); + + Intrepid2::CellTools cell_tools; + MDFieldArrayFactory af(prefix_,true); + + const shards::CellTopology & cell_topology = *(int_rule->topology); + const int cell_dim = cell_topology.getDimension(); + const int subcell_dim = cell_topology.getDimension()-1; + const int num_subcells = cell_topology.getSubcellCount(subcell_dim); + const int num_space_dim = int_rule->topology->getDimension(); + const int num_ip = int_rule->num_points; + const int num_points_on_face = num_ip / num_subcells; + + auto aux = af.template buildStaticArray("surface_normals",num_cells_,num_ip,num_space_dim); + + // We only need the jacobian if we're not in 1D + ConstArray_CellIPDimDim jacobian; + if(cell_dim != 1) + jacobian = getJacobian(false,force); + + // We get to build up our surface normals one 'side' at a time + for(int subcell_index=0; subcell_index("side_normals",num_evaluate_cells_,num_points_on_face,cell_dim); + if(cell_dim == 1){ + + const int other_subcell_index = (subcell_index==0) ? 1 : 0; + auto in_node_coordinates_k = getNodeCoordinates().get_view(); + const auto min = std::numeric_limits::type>::min(); + + Kokkos::parallel_for("compute normals 1D",num_evaluate_cells_,KOKKOS_LAMBDA (const int cell) { + Scalar norm = (in_node_coordinates_k(cell,subcell_index,0) - in_node_coordinates_k(cell,other_subcell_index,0)); + side_normals(cell,0,0) = norm / fabs(norm + min); + }); + + } else { + + // Iterating over side points + Kokkos::MDRangePolicy> policy({0,0},{num_evaluate_cells_,num_points_on_face}); + + auto side_jacobian = Kokkos::DynRankView("side_jac",num_evaluate_cells_,num_points_on_face,cell_dim,cell_dim); + Kokkos::parallel_for("Copy jacobian to side jacobian",policy,KOKKOS_LAMBDA (const int cell,const int point) { + for(int dim=0;dim 0.){ + n = Kokkos::Experimental::sqrt(n); + for(int dim=0;dim> policy({0,0},{num_evaluate_cells_,num_points_on_face}); + Kokkos::parallel_for("copy surface normals", policy,KOKKOS_LAMBDA (const int cell,const int point) { + for(int dim=0;dim 0) + correctVirtualNormals(aux, num_evaluate_cells_ - num_virtual_cells_, num_virtual_cells_, *int_rule->topology, *side_connectivity_); + } + + if(cache){ + surface_normals = aux; + surface_normals_evaluated_ = true; + } + + return aux; + +} + +template +typename IntegrationValues2::ConstArray_CellIPDimDim +IntegrationValues2:: +getSurfaceRotationMatrices(const bool cache, + const bool force) const +{ + + if(surface_rotation_matrices_evaluated_ and not force) + return surface_rotation_matrices; + + MDFieldArrayFactory af(prefix_,true); + + const int num_ip = int_rule->num_points; + const int cell_dim = int_rule->topology->getDimension(); + + // This call will handle all the error checking + auto normals = getSurfaceNormals(false,force).get_static_view(); + auto aux = af.template buildStaticArray("surface_rotation_matrices",num_cells_, num_ip, 3, 3); + + Kokkos::MDRangePolicy> policy({0,0},{num_evaluate_cells_,num_ip}); + Kokkos::parallel_for("create surface rotation matrices",policy,KOKKOS_LAMBDA (const int cell,const int point) { + Scalar normal[3]; + for(int i=0;i<3;i++) + normal[i]=0.; + for(int dim=0; dim 0) + correctVirtualRotationMatrices(aux, num_evaluate_cells_ - num_virtual_cells_, num_virtual_cells_, *int_rule->topology, *side_connectivity_); + } + + if(cache){ + surface_rotation_matrices = aux; + surface_rotation_matrices_evaluated_ = true; + } + + return aux; +} + +template +typename IntegrationValues2::ConstArray_CellIPDimDim +IntegrationValues2:: +getCovarientMatrix(const bool cache, + const bool force) const +{ + + if(covarient_evaluated_ and not force) + return covarient; + + MDFieldArrayFactory af(prefix_,true); + + const int num_space_dim = int_rule->topology->getDimension(); + const int num_ip = int_rule->num_points; + + auto jacobian = getJacobian(false,force).get_static_view(); + auto aux = af.template buildStaticArray("covarient",num_cells_, num_ip, num_space_dim,num_space_dim); + + Kokkos::MDRangePolicy> policy({0,0},{num_evaluate_cells_,num_ip}); + Kokkos::parallel_for("evalaute covarient metric tensor",policy,KOKKOS_LAMBDA (const int cell,const int ip) { + // g^{ij} = \frac{\parital x_i}{\partial \chi_\alpha}\frac{\parital x_j}{\partial \chi_\alpha} + for (int i = 0; i < num_space_dim; ++i) { + for (int j = 0; j < num_space_dim; ++j) { + Scalar value(0.0); + for (int alpha = 0; alpha < num_space_dim; ++alpha) + value += jacobian(cell,ip,i,alpha) * jacobian(cell,ip,j,alpha); + aux(cell,ip,i,j) = value; + } + } + }); + PHX::Device::execution_space().fence(); + + if(cache){ + covarient = aux; + covarient_evaluated_ = true; + } + + return aux; + +} + +template +typename IntegrationValues2::ConstArray_CellIPDimDim +IntegrationValues2:: +getContravarientMatrix(const bool cache, + const bool force) const +{ + + if(contravarient_evaluated_ and not force) + return contravarient; + + MDFieldArrayFactory af(prefix_,true); + + const int num_space_dim = int_rule->topology->getDimension(); + const int num_ip = int_rule->num_points; + + auto cov = getCovarientMatrix(false,force); + auto aux = af.template buildStaticArray("contravarient",num_cells_, num_ip, num_space_dim,num_space_dim); + + const auto cell_range = std::make_pair(0,num_evaluate_cells_); + auto s_contravarient = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL,Kokkos::ALL,Kokkos::ALL); + auto s_covarient = Kokkos::subview(cov.get_view(), cell_range,Kokkos::ALL,Kokkos::ALL,Kokkos::ALL); + + Intrepid2::RealSpaceTools::inverse(s_contravarient, s_covarient); + PHX::Device::execution_space().fence(); + + if(cache){ + contravarient = aux; + contravarient_evaluated_ = true; + } + + return aux; + +} + +template +typename IntegrationValues2::ConstArray_CellIP +IntegrationValues2:: +getNormContravarientMatrix(const bool cache, + const bool force) const +{ + + if(norm_contravarient_evaluated_ and not force) + return norm_contravarient; + + MDFieldArrayFactory af(prefix_,true); + + const int num_space_dim = int_rule->topology->getDimension(); + const int num_ip = int_rule->num_points; + + auto con = getContravarientMatrix(false,force).get_static_view(); + auto aux = af.template buildStaticArray("norm_contravarient",num_cells_, num_ip); + + // norm of g_ij + Kokkos::MDRangePolicy> policy({0,0},{num_evaluate_cells_,num_ip}); + Kokkos::parallel_for("evaluate norm_contravarient",policy,KOKKOS_LAMBDA (const int cell,const int ip) { + aux(cell,ip) = 0.0; + for (int i = 0; i < num_space_dim; ++i) { + for (int j = 0; j < num_space_dim; ++j) { + aux(cell,ip) += con(cell,ip,i,j) * con(cell,ip,i,j); + } + } + aux(cell,ip) = Kokkos::Experimental::sqrt(aux(cell,ip)); + }); + PHX::Device::execution_space().fence(); + + if(cache){ + norm_contravarient = aux; + norm_contravarient_evaluated_ = true; + } + + return aux; + +} + +template +typename IntegrationValues2::ConstArray_CellIPDim +IntegrationValues2:: +getCubaturePoints(const bool cache, + const bool force) const +{ + + if(ip_coordinates_evaluated_ and not force) + return ip_coordinates; + + MDFieldArrayFactory af(prefix_,true); + + const int num_space_dim = int_rule->topology->getDimension(); + const int num_ip = int_rule->num_points; + + auto aux = af.template buildStaticArray("ip_coordinates",num_cells_, num_ip, num_space_dim); + + using ID=panzer::IntegrationDescriptor; + const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY); + + auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates()); + + if(is_cv){ + + // CV integration uses a single call to map from physical space to the weighted measure - I assume this is slower than what we do with non-cv integration methods + const auto cell_range = std::make_pair(0,num_evaluate_cells_); + auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL()); + auto s_cub_points = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL()); + + // TODO: We need to pull this apart for control volumes. Right now we calculate both weighted measures/norms and cubature points at the same time + if(int_rule->cv_type == "side"){ + auto scratch = af.template buildStaticArray("scratch",num_evaluate_cells_,num_ip,num_space_dim); + intrepid_cubature->getCubature(s_cub_points,scratch.get_view(),s_node_coord); + } else { + // I think boundary is included as a weighted measure because it has a side embedded in intrepid_cubature + TEUCHOS_ASSERT((int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_BOUNDARY)); + auto scratch = af.template buildStaticArray("scratch",num_evaluate_cells_,num_ip); + intrepid_cubature->getCubature(s_cub_points,scratch.get_view(),s_node_coord); + } + + } else { + + // Don't forget that since we are not caching this, we have to make sure the managed view remains alive while we use the non-const wrapper + auto const_ref_coord = getCubaturePointsRef(false,force); + auto ref_coord = PHX::getNonConstDynRankViewFromConstMDField(const_ref_coord); + + const auto cell_range = std::make_pair(0,num_evaluate_cells_); + auto s_ref_coord = Kokkos::subview(ref_coord, cell_range,Kokkos::ALL(),Kokkos::ALL()); + auto s_coord = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL()); + auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL()); + + Intrepid2::CellTools cell_tools; + cell_tools.mapToPhysicalFrame(s_coord, s_ref_coord, s_node_coord, *(int_rule->topology)); + + } + + PHX::Device::execution_space().fence(); + + if(cache){ + ip_coordinates = aux; + ip_coordinates_evaluated_ = true; + } + + return aux; + +} + + +template +typename IntegrationValues2::ConstArray_CellIPDim +IntegrationValues2:: +getCubaturePointsRef(const bool cache, + const bool force) const +{ + + if(ref_ip_coordinates_evaluated_) + return ref_ip_coordinates; + + using ID=panzer::IntegrationDescriptor; + const bool is_surface = int_rule->getType() == ID::SURFACE; + const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY); + + const int num_space_dim = int_rule->topology->getDimension(); + const int num_ip = int_rule->num_points; + + MDFieldArrayFactory af(prefix_,true); + + Intrepid2::CellTools cell_tools; + + auto aux = af.template buildStaticArray("ref_ip_coordinates",num_cells_, num_ip, num_space_dim); + + if(is_cv){ + + // Control volume reference points are actually generated from the physical points (i.e. reverse to everything else) + + auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates()); + + // Don't forget that since we are not caching this, we have to make sure the managed view remains alive while we use the non-const wrapper + auto const_coord = getCubaturePoints(false,force); + auto coord = PHX::getNonConstDynRankViewFromConstMDField(const_coord); + + const auto cell_range = std::make_pair(0,num_evaluate_cells_); + auto s_ref_coord = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL()); + auto s_coord = Kokkos::subview(coord, cell_range,Kokkos::ALL(),Kokkos::ALL()); + auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL()); + + cell_tools.mapToReferenceFrame(s_ref_coord, s_coord, s_node_coord, *(int_rule->topology)); + + } else if(is_surface){ + + const auto & cell_topology = *int_rule->topology; + const int cell_dim = cell_topology.getDimension(); + const int num_sides = (cell_dim==1) ? 2 : cell_topology.getSideCount(); + const int subcell_dim = cell_dim-1; + const int num_points_on_face = num_ip / num_sides; + const int order = int_rule->getOrder(); + + // Scratch space for storing the points for each side of the cell + auto side_cub_points = af.template buildStaticArray("side_cub_points",num_points_on_face,cell_dim); + + Intrepid2::DefaultCubatureFactory cubature_factory; + + // We get to build up our cubature one side at a time + for(int side=0; side(face_topology,order); + auto tmp_side_cub_weights = Kokkos::DynRankView("tmp_side_cub_weights",num_points_on_face); + auto tmp_side_cub_points = Kokkos::DynRankView("tmp_side_cub_points",num_points_on_face,subcell_dim); + + // Get the reference face points + ic->getCubature(tmp_side_cub_points, tmp_side_cub_weights); + + // Convert from reference face points to reference cell points + cell_tools.mapToReferenceSubcell(side_cub_points.get_view(), tmp_side_cub_points, subcell_dim, side, cell_topology); + } + + PHX::Device::execution_space().fence(); + + // Copy from the side allocation to the surface allocation + Kokkos::MDRangePolicy> policy({0,0,0},{num_evaluate_cells_,num_points_on_face, num_space_dim}); + Kokkos::parallel_for("copy values",policy,KOKKOS_LAMBDA (const int cell,const int point, const int dim) { + aux(cell,point_offset + point,dim) = side_cub_points(point,dim); + }); + } + + } else { + + // Haven't set this up yet + TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->isSide() && num_space_dim==1, + "ERROR: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do " + << "non-natural integration rules."); + + auto cub_points = getUniformCubaturePointsRef(false,force,false); + + Kokkos::MDRangePolicy> policy({0,0,0},{num_evaluate_cells_,num_ip,num_space_dim}); + Kokkos::parallel_for(policy, KOKKOS_LAMBDA(const int & cell, const int & ip, const int & dim){ + aux(cell,ip,dim) = cub_points(ip,dim); + }); + } + + PHX::Device::execution_space().fence(); + + if(requires_permutation_) + applyPermutation(aux, permutations_); + + if(cache){ + ref_ip_coordinates = aux; + ref_ip_coordinates_evaluated_ = true; + } + + return aux; + +} + +template +void +IntegrationValues2:: +evaluateEverything() +{ + + using ID=panzer::IntegrationDescriptor; + const bool is_surface = int_rule->getType() == ID::SURFACE; + const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY); + const bool is_side = int_rule->isSide(); + + resetArrays(); + + // Base cubature stuff + if(is_cv){ + getCubaturePoints(true,true); + getCubaturePointsRef(true,true); + } else { + if(not is_surface){ + getUniformCubaturePointsRef(true,true); + getUniformCubatureWeightsRef(true,true); + if(is_side) + getUniformSideCubaturePointsRef(true,true); + } + getCubaturePointsRef(true,true); + getCubaturePoints(true,true); + } + + // Measure stuff + getJacobian(true,true); + getJacobianDeterminant(true,true); + getJacobianInverse(true,true); + if(int_rule->cv_type == "side") + getWeightedNormals(true,true); + else + getWeightedMeasure(true,true); + + // Surface stuff + if(is_surface){ + getSurfaceNormals(true,true); + getSurfaceRotationMatrices(true,true); + } + + // Stabilization stuff + if(not (is_surface or is_cv)){ + getContravarientMatrix(true,true); + getCovarientMatrix(true,true); + getNormContravarientMatrix(true,true); + } - cell_tools.setJacobianDet(s_jac_det, s_jac); } #define INTEGRATION_VALUES2_INSTANTIATION(SCALAR) \ diff --git a/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.hpp b/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.hpp index bab023aa6e54..4b0d9128ed3d 100644 --- a/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.hpp +++ b/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.hpp @@ -57,31 +57,6 @@ namespace panzer { class SubcellConnectivity; - /** - * \brief Swap the ordering of quadrature points in a specified cell. - * - * \param[in] cell Cell index - * \param[in] a Quadrature point a - * \param[in] b Quadrature point b - * - * NOTE: this should be a member function of IntegrationValues but - * lambda capture generates a ton of cuda compiler warnings. Making - * it a stand alone function fixes this. - */ - template - KOKKOS_INLINE_FUNCTION - void swapQuadraturePoints(int cell,int a,int b, - T0& ref_ip_coordinates, - T1& ip_coordinates, - T2& weighted_measure, - T3& jac, - T4& jac_det, - T5& jac_inv, - T6& surface_normals, - T7& surface_rotation_matrices); - template class IntegrationValues2 { public: @@ -90,24 +65,38 @@ namespace panzer { typedef PHX::MDField ArrayDynamic; typedef PHX::MDField DblArrayDynamic; - typedef PHX::MDField Array_IP; - typedef PHX::MDField Array_IPDim; - - typedef PHX::MDField Array_Point; - typedef PHX::MDField Array_CellIP; - typedef PHX::MDField Array_CellIPDim; - typedef PHX::MDField Array_CellIPDimDim; + typedef PHX::MDField Array_IP; + typedef PHX::MDField Array_IPDim; + typedef PHX::MDField Array_Point; + typedef PHX::MDField Array_CellIP; + typedef PHX::MDField Array_CellIPDim; + typedef PHX::MDField Array_CellIPDimDim; + typedef PHX::MDField Array_CellBASISDim; + + typedef PHX::MDField ConstArray_IP; + typedef PHX::MDField ConstArray_IPDim; + typedef PHX::MDField ConstArray_Point; + typedef PHX::MDField ConstArray_CellIP; + typedef PHX::MDField ConstArray_CellIPDim; + typedef PHX::MDField ConstArray_CellIPDimDim; + typedef PHX::MDField ConstArray_CellBASISDim; + + /** + * \brief Base constructor + * + * \param[in] pre Prefix to apply to all internal field names + * \param[in] allocArrays (Classic Interface Only) Allocate array data in 'setupArrays' call + */ + IntegrationValues2(const std::string & pre="", + const bool allocArrays=false); - typedef PHX::MDField Array_CellBASISDim; - IntegrationValues2(const std::string & pre="",bool allocArrays=false) - : alloc_arrays(allocArrays), prefix(pre), ddims_(1,0) {} + // ===================================================================================================== + // Classic Interface (DEPRECATED) //! Sizes/allocates memory for arrays void setupArrays(const Teuchos::RCP& ir); - void setupArraysForNodeRule(const Teuchos::RCP& ir); - /** \brief Evaluate basis values. @param vertex_coordinates [in] Cell vertex coordinates, not @@ -121,7 +110,8 @@ namespace panzer { */ void evaluateValues(const PHX::MDField & vertex_coordinates, const int num_cells = -1, - const Teuchos::RCP & face_connectivity = Teuchos::null); + const Teuchos::RCP & face_connectivity = Teuchos::null, + const int num_virtual_cells = -1); /** \brief Match IP. @@ -141,60 +131,382 @@ namespace panzer { const PHX::MDField & other_ip_coordinates, const int num_cells = -1); - Array_IPDim cub_points; // - Array_IPDim side_cub_points; // points on face topology (dim-1) - Array_IP cub_weights; // - Array_CellBASISDim node_coordinates; // - Array_CellIPDimDim jac; // - Array_CellIPDimDim jac_inv; // - Array_CellIP jac_det; // - Array_CellIP weighted_measure; // - Array_CellIPDim weighted_normals; // - - Array_CellIPDim surface_normals; // - Array_CellIPDimDim surface_rotation_matrices; // + // Reference space quantities + mutable Array_IPDim cub_points; // + mutable Array_IPDim side_cub_points; // points on face topology (dim-1) + mutable Array_IP cub_weights; // + + // Physical space quantities + mutable Array_CellBASISDim node_coordinates; // + mutable Array_CellIPDimDim jac; // + mutable Array_CellIPDimDim jac_inv; // + mutable Array_CellIP jac_det; // + mutable Array_CellIP weighted_measure; // + mutable Array_CellIPDim weighted_normals; // + mutable Array_CellIPDim surface_normals; // + mutable Array_CellIPDimDim surface_rotation_matrices; // // this (appears) is a matrix where the first row is the "normal" direction // and the remaining two rows lie in the hyperplane - Teuchos::RCP int_rule; - - Teuchos::RCP> intrepid_cubature; - // for Shakib stabilization - Array_CellIPDimDim covarient; - Array_CellIPDimDim contravarient; - Array_CellIP norm_contravarient; + mutable Array_CellIPDimDim covarient; + mutable Array_CellIPDimDim contravarient; + mutable Array_CellIP norm_contravarient; // integration points - Array_CellIPDim ip_coordinates; // - Array_CellIPDim ref_ip_coordinates; // for Control Volumes or Surface integrals + mutable Array_CellIPDim ip_coordinates; // + mutable Array_CellIPDim ref_ip_coordinates; // for Control Volumes or Surface integrals - DblArrayDynamic dyn_cub_points, dyn_side_cub_points, dyn_cub_weights; - DblArrayDynamic dyn_phys_cub_points, dyn_phys_cub_weights, dyn_phys_cub_norms, dyn_node_coordinates; - Array_Point scratch_for_compute_side_measure; // size: span() == jac.span() + Teuchos::RCP int_rule; + + Teuchos::RCP> intrepid_cubature; + + // ===================================================================================================== + // Lazy evaluation interface + + /** + * The lazy evaluation construction path is designed such that you only allocate and fill arrays on demand, + * with an option of caching those generated fields. This is useful for when we are worried about a + * code's memory footprint. + * + * The process for setting up one of these objects is to initialize the IntegrationValues2 class as follows: + * + * IntegrationValues2 values; // Constructor + * values.setup(ir, node_coordinates); // Required - must come before all other calls + * values.setupPermutations(...); // Optional - only needed for surface integration and some side integration rules - must come before get*() calls + * auto array = values.get*(); // Lazy evaluate whatever field you need + */ + + /** + * \brief Main setup call for the lazy evaluation interface + * + * \todo Instead of IntegrationRule, we just need to load the integration descriptor and the cell topology + * + * \param[in] ir Integration rule descripting integration scheme + * \param[in] node_coordinates Node/Vertex coordinates describing cell geometry + * \param[in] num_cells In case you need to only generate integration values for the first 'num_cells' of the node_coordinates - defaults to all cells + */ + void + setup(const Teuchos::RCP& ir, + const PHX::MDField & node_coordinates, + const int num_cells = -1); + + /** + * \brief Initialize the permutation arrays given a face connectivity + * + * \note REQUIRED FOR SURFACE INTEGRATION + * \note Must be called AFTER setup + * \note Virtual cells have a unique way to generate surface normals and rotation matrices, hence we need to know how many are included + * + * \param[in] face_connectivity Connectivity describing how sides are connected to cells + * \param[in] num_virtual_cells Number of virtual cells included in the node coordinates (found at end of node coordinate array) + */ + void + setupPermutations(const Teuchos::RCP & face_connectivity, + const int num_virtual_cells); + + /** + * \brief Initialize the permutation arrays given another IntegrationValues2::getCubaturePoints() array + * + * \note Required if you want points to line up between two integration values (e.g. for side integration) + * \note Must be called AFTER setup + * + * \param[in] other_ip_coordinates Cubature points to align with + */ + void + setupPermutations(const PHX::MDField & other_ip_coordinates); + + /** + * \brief Get the uniform cubature points + * + * \note DEPRECATED Please use getCubaturePointsRef call instead + * \note Option is only supported for volume integration, and even then, may not align with the getCubaturePointsRef call + * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * \param[in] apply_permutation ADVANCED Do not change this unless you know what it does (it can break things) + * + * \return Array + */ + ConstArray_IPDim + getUniformCubaturePointsRef(const bool cache = true, + const bool force = false, + const bool apply_permutation = true) const; + + /** + * \brief Get the uniform cubature points for a side + * + * \note DEPRECATED Please use getCubaturePointsRef call instead + * \note Option is only supported for side integration, and even then, may not align with the getCubaturePointsRef call + * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * \param[in] apply_permutation ADVANCED Do not change this unless you know what it does (it can break things) + * + * \return Array + */ + ConstArray_IPDim + getUniformSideCubaturePointsRef(const bool cache = true, + const bool force = false, + const bool apply_permutation = true) const; + + /** + * \brief Get the uniform cubature weights + * + * \note DEPRECATED Please do not use + * \note Option is only supported for some integration types + * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * \param[in] apply_permutation ADVANCED Do not change this unless you know what it does (it can break things) + * + * \return Array + */ + ConstArray_IP + getUniformCubatureWeightsRef(const bool cache = true, + const bool force = false, + const bool apply_permutation = true) const; - /// This should be a private method, but using lambdas on cuda forces this to be public. - void evaluateRemainingValues(const PHX::MDField & in_node_coordinates, const int in_num_cells); - /// This should be a private method, but using lambdas on cuda forces this to be public. - void getCubatureCV(const PHX::MDField & in_node_coordinates, const int in_num_cells); + /** + * \brief Get the node coordinates describing the geometry of the mesh + * + * \return Array + */ + ConstArray_CellBASISDim + getNodeCoordinates() const; + + /** + * \brief Get the Jacobian matrix evaluated at the cubature points + * + * \note Support: VOLUME, SURFACE, SIDE, CV_VOLUME, CV_SIDE, CV_BOUNDARY + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIPDimDim + getJacobian(const bool cache = true, + const bool force = false) const; + + + /** + * \brief Get the inverse of the Jacobian matrix evaluated at the cubature points + * + * \note Support: VOLUME, SURFACE, SIDE, CV_VOLUME, CV_SIDE, CV_BOUNDARY + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIPDimDim + getJacobianInverse(const bool cache = true, + const bool force = false) const; + + /** + * \brief Get the determinant of the Jacobian matrix evaluated at the cubature points + * + * \note Support: VOLUME, SURFACE, SIDE, CV_VOLUME, CV_SIDE, CV_BOUNDARY + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIP + getJacobianDeterminant(const bool cache = true, + const bool force = false) const; + + /** + * \brief Get the weighted measure (integration weights) + * + * \note Support: VOLUME, SURFACE, SIDE, CV_VOLUME, CV_BOUNDARY + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIP + getWeightedMeasure(const bool cache = true, + const bool force = false) const; + + /** + * \brief Get the weighted normals + * + * \note Support: CV_SIDE + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIPDim + getWeightedNormals(const bool cache = true, + const bool force = false) const; + + /** + * \brief Get the surface normals + * + * \note Support: SURFACE + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIPDim + getSurfaceNormals(const bool cache = true, + const bool force = false) const; + + /** + * \brief Get the surface rotation matrices + * + * \note Support: SURFACE + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIPDimDim + getSurfaceRotationMatrices(const bool cache = true, + const bool force = false) const; + + /** + * \brief Get the covarient matrix + * + * cov(i,j) = jacobian(i,k) * jacobian(j,k) + * + * \note Support: VOLUME, SURFACE, SIDE, CV_VOLUME, CV_SIDE, CV_BOUNDARY + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIPDimDim + getCovarientMatrix(const bool cache = true, + const bool force = false) const; + + /** + * \brief Get the contravarient matrix + * + * contra = (getCovarientMatrix())^{-1} + * + * \note Support: VOLUME, SURFACE, SIDE, CV_VOLUME, CV_SIDE, CV_BOUNDARY + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIPDimDim + getContravarientMatrix(const bool cache = true, + const bool force = false) const; + + /** + * \brief Get the contravarient matrix + * + * norm = sqrt(\sum_{ij} cov(i,j) * cov(i,j)) + * + * \note Support: VOLUME, SURFACE, SIDE, CV_VOLUME, CV_SIDE, CV_BOUNDARY + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIP + getNormContravarientMatrix(const bool cache = true, + const bool force = false) const; + + /** + * \brief Get the cubature points in physical space + * + * \note Support: VOLUME, SURFACE, SIDE, CV_VOLUME, CV_SIDE, CV_BOUNDARY + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIPDim + getCubaturePoints(const bool cache = true, + const bool force = false) const; + + /** + * \brief Get the cubature points in the reference space + * + * \note Support: VOLUME, SURFACE, SIDE, CV_VOLUME, CV_SIDE, CV_BOUNDARY + * * + * \param[in] cache If true, the result will be stored in the IntegrationValues2 class + * \param[in] force Force the re-evaluation of the array + * + * \return Array + */ + ConstArray_CellIPDim + getCubaturePointsRef(const bool cache = true, + const bool force = false) const; - /// This should be a private method, but using lambdas on cuda forces this to be public. - void generateSurfaceCubatureValues(const PHX::MDField & in_node_coordinates, const int in_num_cells,const SubcellConnectivity & face_connectivity); + // ===================================================================================================== protected: - // TODO: Make this a utility function that only exists in source file - Teuchos::RCP> getIntrepidCubature(const panzer::IntegrationRule & ir) const; + // Reset all the lazy evaluation arrays + void + resetArrays(); + + // Number of cells in mesh + int num_cells_; + + // Number of cells in mesh to evaluate + int num_evaluate_cells_; + + // Number of virtual cells in the mesh - used for surface evaluations + int num_virtual_cells_; + + // Permutations (used to re-orient arrays similar to orientations in BasisValues2) + bool requires_permutation_; + + // Array contains the mapping from uniform reference space to permuted space + PHX::MDField permutations_; + + // TODO: There is a way around this, but it will require some work + // Subcell connectivity is required for surface evaluations (normals and rotation matrices) + Teuchos::RCP side_connectivity_; + + // Lazy evaluation checks + mutable bool cub_points_evaluated_; + mutable bool side_cub_points_evaluated_; + mutable bool cub_weights_evaluated_; + mutable bool node_coordinates_evaluated_; + mutable bool jac_evaluated_; + mutable bool jac_inv_evaluated_; + mutable bool jac_det_evaluated_; + mutable bool weighted_measure_evaluated_; + mutable bool weighted_normals_evaluated_; + mutable bool surface_normals_evaluated_; + mutable bool surface_rotation_matrices_evaluated_; + mutable bool covarient_evaluated_; + mutable bool contravarient_evaluated_; + mutable bool norm_contravarient_evaluated_; + mutable bool ip_coordinates_evaluated_; + mutable bool ref_ip_coordinates_evaluated_; + + // Backward compatibility call that evaluates all internal values for CV, surface, side, or volume integration schemes + void + evaluateEverything(); private: - bool alloc_arrays; - std::string prefix; + + bool alloc_arrays_; + std::string prefix_; std::vector ddims_; - void getCubature(const PHX::MDField & in_node_coordinates, const int in_num_cells); - void evaluateValuesCV(const PHX::MDField & vertex_coordinates,const int in_num_cells); }; } // namespace panzer diff --git a/packages/panzer/disc-fe/src/Panzer_IntegrationValues2_impl.hpp b/packages/panzer/disc-fe/src/Panzer_IntegrationValues2_impl.hpp deleted file mode 100644 index 2af2de2e8381..000000000000 --- a/packages/panzer/disc-fe/src/Panzer_IntegrationValues2_impl.hpp +++ /dev/null @@ -1,93 +0,0 @@ -#ifndef PANZER_INTEGRATION_VALUES_2_IMPL_HPP -#define PANZER_INTEGRATION_VALUES_2_IMPL_HPP - -/* Normally, we would fold this definition into the .cpp file since we - use ETI to compile this object. However, for unit testing the stand - alone device function swapQuadraturePoints, we need the definition - at the calling point to avoid requiring relocatable device code for - cuda. So for the IntegrationValues2 object, if there is a - standalone device member function, put it in this file (for unit - testing). - */ -namespace panzer { - -template -void -swapQuadraturePoints(int cell, - int a, - int b, - T0& ref_ip_coordinates, - T1& ip_coordinates, - T2& weighted_measure, - T3& jac, - T4& jac_det, - T5& jac_inv, - T6& surface_normals, - T7& surface_rotation_matrices) -{ - const int new_cell_point = a; - const int old_cell_point = b; - - const int cell_dim = ref_ip_coordinates.extent(2); - -#ifdef PANZER_DEBUG - KOKKOS_ASSERT(cell < ip_coordinates.extent_int(0)); - KOKKOS_ASSERT(a < ip_coordinates.extent_int(1)); - KOKKOS_ASSERT(b < ip_coordinates.extent_int(1)); - KOKKOS_ASSERT(cell >= 0); - KOKKOS_ASSERT(a >= 0); - KOKKOS_ASSERT(b >= 0); -#endif - - // If this is a DFAD type, we will need to fix allocation to size the derivative array. - Scalar hold; - - hold = weighted_measure(cell,new_cell_point); - weighted_measure(cell,new_cell_point) = weighted_measure(cell,old_cell_point); - weighted_measure(cell,old_cell_point) = hold; - - hold = jac_det(cell,new_cell_point); - jac_det(cell,new_cell_point) = jac_det(cell,old_cell_point); - jac_det(cell,old_cell_point) = hold; - - for(int dim=0;dim > iv, - const int num_real_cells, - const int num_virtual_cells, - const Teuchos::RCP cell_topology, - const SubcellConnectivity & face_connectivity) -{ - - if (iv->int_rule->getType() != panzer::IntegrationDescriptor::SURFACE) - return; - - // IntegrationValues2 doesn't know anything about virtual cells, so it sets up incorrect normals for those. - // What we want is for the adjoining face of the virtual cell to have normals that are the negated real cell's normals. - // we correct the normals here: - const int space_dim = cell_topology->getDimension(); - const int faces_per_cell = cell_topology->getSubcellCount(space_dim-1); - const int points = iv->surface_normals.extent_int(1); - const int points_per_face = points / faces_per_cell; - - auto surface_normals_view = PHX::as_view(iv->surface_normals); - auto surface_normals_h = Kokkos::create_mirror_view(surface_normals_view); - Kokkos::deep_copy(surface_normals_h, surface_normals_view); - - auto surface_rotation_matrices_view = PHX::as_view(iv->surface_rotation_matrices); - auto surface_rotation_matrices_h = Kokkos::create_mirror_view(surface_rotation_matrices_view); - Kokkos::deep_copy(surface_rotation_matrices_h, surface_rotation_matrices_view); - - - for (int virtual_cell_ordinal=0; virtual_cell_ordinal= 0) - { - virtual_local_face_id = local_face_id; - break; - } - } - if (face_ordinal >= 0) - { - const int first_cell_for_face = face_connectivity.cellForSubcellHost(face_ordinal, 0); - const panzer::LocalOrdinal other_side = (first_cell_for_face == virtual_cell) ? 1 : 0; - const panzer::LocalOrdinal real_cell = face_connectivity.cellForSubcellHost(face_ordinal,other_side); - const panzer::LocalOrdinal face_in_real_cell = face_connectivity.localSubcellForSubcellHost(face_ordinal,other_side); - TEUCHOS_ASSERT(real_cell < num_real_cells); - for (int point_ordinal=0; point_ordinal & WorksetDetails:: -getIntegrationValues(const panzer::IntegrationDescriptor & description) const +getIntegrationValues(const panzer::IntegrationDescriptor & description, + const bool lazy_version) const { TEUCHOS_ASSERT(setup_); + // We need unique keys for the lazy copy or else we get some weird behavior + size_t key = description.getKey(); + if(lazy_version) + panzer::hash_combine(key, 123); + // Check if exists - const auto itr = integration_values_map_.find(description.getKey()); + const auto itr = integration_values_map_.find(key); if(itr != integration_values_map_.end()) return *(itr->second); @@ -345,21 +247,33 @@ getIntegrationValues(const panzer::IntegrationDescriptor & description) const TEUCHOS_TEST_FOR_EXCEPT_MSG(description.getSide() != getSubcellIndex(), "Workset::getIntegrationValues : Attempted to build integration values for side '"<("",true)); - iv->setupArrays(ir); - iv->evaluateValues(getCellVertices(), numCells(), face_connectivity_); + // Create the integration values object + Teuchos::RCP> iv; + if(lazy_version){ + iv = Teuchos::rcp(new IntegrationValues2()); + + iv->setup(ir,getCellVertices(),numCells()); + + // Surface integration schemes need to properly "permute" their entries to line up the surface points between cells + if(description.getType() == panzer::IntegrationDescriptor::SURFACE) + iv->setupPermutations(face_connectivity_, numVirtualCells()); - correctVirtualNormals(iv, num_owned_cells_ + num_ghost_cells_, num_virtual_cells_, cell_topology_, *face_connectivity_); + } else { + + iv = Teuchos::rcp(new IntegrationValues2("",true)); + iv->setupArrays(ir); + iv->evaluateValues(getCellVertices(), numCells(), face_connectivity_, numVirtualCells()); - // This is an advanced feature that requires changes to the workset construction - // Basically there needs to be a way to grab the side assembly for both "details" belonging to the same workset, which requires a refactor - // Alternatively, we can do this using a face connectivity object, but the refactor is more important atm. - TEUCHOS_ASSERT(not (options_.side_assembly_ and options_.align_side_points_)); + // This is an advanced feature that requires changes to the workset construction + // Basically there needs to be a way to grab the side assembly for both "details" belonging to the same workset, which requires a refactor + // Alternatively, we can do this using a face connectivity object, but the refactor is more important atm. + TEUCHOS_ASSERT(not (options_.side_assembly_ and options_.align_side_points_)); + + } - integration_values_map_[description.getKey()] = iv; + integration_values_map_[key] = iv; ir_degrees->push_back(iv->int_rule->cubature_degree); int_rules.push_back(iv); @@ -406,18 +320,28 @@ getBasisValues(const panzer::BasisDescriptor & basis_description, { TEUCHOS_ASSERT(setup_); + // We need unique keys for the lazy copy or else we get some weird behavior + size_t basis_key = basis_description.getKey(); + if(lazy_version) + panzer::hash_combine(basis_key, 123); + + // We need unique keys for the lazy copy or else we get some weird behavior + size_t integration_key = integration_description.getKey(); + if(lazy_version) + panzer::hash_combine(integration_key, 123); + // Check if exists - const auto itr = basis_integration_values_map_.find(basis_description.getKey()); + const auto itr = basis_integration_values_map_.find(basis_key); if(itr != basis_integration_values_map_.end()){ const auto & submap = itr->second; - const auto itr2 = submap.find(integration_description.getKey()); + const auto itr2 = submap.find(integration_key); if(itr2 != submap.end()) return *(itr2->second); } // Get the integration values for this description - const auto & iv = getIntegrationValues(integration_description); + const auto & iv = getIntegrationValues(integration_description,lazy_version); auto bir = Teuchos::rcp(new BasisIRLayout(basis_description.getType(), basis_description.getOrder(), *iv.int_rule)); Teuchos::RCP> biv; @@ -429,12 +353,12 @@ getBasisValues(const panzer::BasisDescriptor & basis_description, biv = Teuchos::rcp(new BasisValues2()); if(integration_description.getType() == IntegrationDescriptor::VOLUME) - biv->setupUniform(bir, iv.cub_points, iv.jac, iv.jac_det, iv.jac_inv); + biv->setupUniform(bir, iv.getUniformCubaturePointsRef(false), iv.getJacobian(false), iv.getJacobianDeterminant(false), iv.getJacobianInverse(false)); else - biv->setup(bir, iv.ref_ip_coordinates, iv.jac, iv.jac_det, iv.jac_inv); + biv->setup(bir, iv.getCubaturePointsRef(false), iv.getJacobian(false), iv.getJacobianDeterminant(false), iv.getJacobianInverse(false)); biv->setOrientations(options_.orientations_, numOwnedCells()+numGhostCells()); - biv->setWeightedMeasure(iv.weighted_measure); + biv->setWeightedMeasure(iv.getWeightedMeasure(false)); biv->setCellVertexCoordinates(cell_vertex_coordinates); } else { @@ -486,7 +410,7 @@ getBasisValues(const panzer::BasisDescriptor & basis_description, } - basis_integration_values_map_[basis_description.getKey()][integration_description.getKey()] = biv; + basis_integration_values_map_[basis_key][integration_key] = biv; bases.push_back(biv); basis_names->push_back(biv->basis_layout->name()); @@ -501,6 +425,7 @@ getPointValues(const panzer::PointDescriptor & description) const { TEUCHOS_ASSERT(setup_); + // Check if exists const auto itr = point_values_map_.find(description.getKey()); if(itr != point_values_map_.end()) @@ -534,8 +459,13 @@ getBasisValues(const panzer::BasisDescriptor & basis_description, { TEUCHOS_ASSERT(setup_); + // We need unique keys for the lazy copy or else we get some weird behavior + size_t basis_key = basis_description.getKey(); + if(lazy_version) + panzer::hash_combine(basis_key, 123); + // Check if exists - const auto itr = basis_point_values_map_.find(basis_description.getKey()); + const auto itr = basis_point_values_map_.find(basis_key); if(itr != basis_point_values_map_.end()){ const auto & submap = itr->second; const auto itr2 = submap.find(point_description.getKey()); @@ -581,7 +511,7 @@ getBasisValues(const panzer::BasisDescriptor & basis_description, } - basis_point_values_map_[basis_description.getKey()][point_description.getKey()] = bpv; + basis_point_values_map_[basis_key][point_description.getKey()] = bpv; return *bpv; diff --git a/packages/panzer/disc-fe/src/Panzer_Workset.hpp b/packages/panzer/disc-fe/src/Panzer_Workset.hpp index df416cb46e81..deb540539135 100644 --- a/packages/panzer/disc-fe/src/Panzer_Workset.hpp +++ b/packages/panzer/disc-fe/src/Panzer_Workset.hpp @@ -236,11 +236,13 @@ namespace panzer { * \throws If setup has not been called * * \param[in] description Descriptor for integration scheme + * \param[in] lazy_version Get an empty IntegrationValues2 object that will construct/allocate itself on demand (less memory - EXPERIMENTAL) * * \return Object containing integration values */ const panzer::IntegrationValues2 & - getIntegrationValues(const panzer::IntegrationDescriptor & description) const; + getIntegrationValues(const panzer::IntegrationDescriptor & description, + const bool lazy_version=false) const; /* * \brief Grab the basis values for a given basis description 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 d08b356c9fd0..8fd5c4aa61a7 100644 --- a/packages/panzer/disc-fe/test/core_tests/basis_values2.cpp +++ b/packages/panzer/disc-fe/test/core_tests/basis_values2.cpp @@ -138,7 +138,7 @@ namespace panzer { int_values.jac, int_values.jac_det, int_values.jac_inv, - int_values.weighted_measure, + int_values.weighted_measure, node_coordinates); out << "check output" << std::endl; diff --git a/packages/panzer/disc-fe/test/core_tests/integration_rule.cpp b/packages/panzer/disc-fe/test/core_tests/integration_rule.cpp index 7e140d520eb0..f1bb0a562780 100644 --- a/packages/panzer/disc-fe/test/core_tests/integration_rule.cpp +++ b/packages/panzer/disc-fe/test/core_tests/integration_rule.cpp @@ -117,7 +117,8 @@ namespace panzer { const int num_cells = 20; const int base_cell_dimension = 2; - const panzer::CellData cell_data(num_cells,topo); + const int cell_local_side_id = 1; + const panzer::CellData cell_data(num_cells,cell_local_side_id,topo); std::string cv_type = "side"; panzer::IntegrationRule int_rule(cell_data,cv_type); diff --git a/packages/panzer/disc-fe/test/core_tests/integration_values2.cpp b/packages/panzer/disc-fe/test/core_tests/integration_values2.cpp index f50ef6961b8f..41d4bd9deef0 100644 --- a/packages/panzer/disc-fe/test/core_tests/integration_values2.cpp +++ b/packages/panzer/disc-fe/test/core_tests/integration_values2.cpp @@ -48,7 +48,6 @@ #include "Panzer_CellData.hpp" #include "Panzer_IntegrationRule.hpp" #include "Panzer_IntegrationValues2.hpp" -#include "Panzer_IntegrationValues2_impl.hpp" // for unit testing wihtout CUDA RDC #include "Panzer_ArrayTraits.hpp" #include "Panzer_CommonArrayFactories.hpp" #include "Panzer_SubcellConnectivity.hpp" @@ -73,11 +72,8 @@ namespace panzer { RCP int_rule = rcp(new IntegrationRule(cubature_degree, cell_data)); - panzer::IntegrationValues2 int_values("prefix_",true); panzer::MDFieldArrayFactory af("prefix_",true); - int_values.setupArrays(int_rule); - const int num_vertices = int_rule->topology->getNodeCount(); PHX::MDField node_coordinates = af.buildStaticArray("nc",num_cells, num_vertices, base_cell_dimension); @@ -107,12 +103,13 @@ namespace panzer { node_coordinates_k(cell,3,y) = 1.0; }); - int_values.evaluateValues(node_coordinates); + panzer::IntegrationValues2 int_values("prefix_"); + int_values.setup(int_rule, node_coordinates); - TEST_EQUALITY(int_values.ip_coordinates.extent(1), 4); + TEST_EQUALITY(int_values.getCubaturePoints().extent(1), 4); double realspace_x_coord = (1.0/std::sqrt(3.0) + 1.0) / 2.0; double realspace_y_coord = (1.0/std::sqrt(3.0) + 1.0) / 2.0; - auto tmp_coords = int_values.ip_coordinates; + auto tmp_coords = int_values.getCubaturePoints(); auto tmp_coords_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),tmp_coords.get_view()); TEST_FLOATING_EQUALITY(tmp_coords_host(0,0,0), realspace_x_coord, 1.0e-8); @@ -121,8 +118,8 @@ namespace panzer { } TEUCHOS_UNIT_TEST(integration_values, control_volume) - { - Teuchos::RCP topo = + { + Teuchos::RCP topo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >())); const int num_cells = 20; @@ -130,24 +127,17 @@ namespace panzer { const panzer::CellData cell_data(num_cells,topo); std::string cv_type = "volume"; - RCP int_rule_vol = + RCP int_rule_vol = rcp(new IntegrationRule(cell_data, cv_type)); - - panzer::IntegrationValues2 int_values_vol("prefix_",true); - panzer::MDFieldArrayFactory af("prefix_",true); - int_values_vol.setupArrays(int_rule_vol); + panzer::MDFieldArrayFactory af("prefix_",true); cv_type = "side"; - RCP int_rule_side = + RCP int_rule_side = rcp(new IntegrationRule(cell_data, cv_type)); - - panzer::IntegrationValues2 int_values_side("prefix_",true); - - int_values_side.setupArrays(int_rule_side); const int num_vertices = int_rule_vol->topology->getNodeCount(); - PHX::MDField node_coordinates + PHX::MDField node_coordinates = af.buildStaticArray("nc",num_cells, num_vertices, base_cell_dimension); // Set up node coordinates. Here we assume the following @@ -175,34 +165,38 @@ namespace panzer { node_coordinates_k(cell,3,y) = 1.0; }); - int_values_vol.evaluateValues(node_coordinates); - int_values_side.evaluateValues(node_coordinates); - - TEST_EQUALITY(int_values_vol.ip_coordinates.extent(1), 4); - TEST_EQUALITY(int_values_side.ip_coordinates.extent(1), 4); - TEST_EQUALITY(int_values_side.weighted_normals.extent(1), 4); + + panzer::IntegrationValues2 int_values_vol("prefix_"); + int_values_vol.setup(int_rule_vol, node_coordinates); + + panzer::IntegrationValues2 int_values_side("prefix_"); + int_values_side.setup(int_rule_side, node_coordinates); + + TEST_EQUALITY(int_values_vol.getCubaturePoints().extent(1), 4); + TEST_EQUALITY(int_values_side.getCubaturePoints().extent(1), 4); + TEST_EQUALITY(int_values_side.getWeightedNormals().extent(1), 4); double realspace_x_coord_1 = 0.25; double realspace_y_coord_1 = 0.25; - auto tmp_coords = int_values_vol.ip_coordinates; + auto tmp_coords = int_values_vol.getCubaturePoints(); auto tmp_coords_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),tmp_coords.get_view()); - TEST_FLOATING_EQUALITY(tmp_coords_host(0,0,0), + TEST_FLOATING_EQUALITY(tmp_coords_host(0,0,0), realspace_x_coord_1, 1.0e-8); - TEST_FLOATING_EQUALITY(tmp_coords_host(0,0,1), + TEST_FLOATING_EQUALITY(tmp_coords_host(0,0,1), realspace_y_coord_1, 1.0e-8); double realspace_x_coord_2 = 0.5; double realspace_y_coord_2 = 0.25; - auto tmp_side_coords = int_values_side.ip_coordinates; + auto tmp_side_coords = int_values_side.getCubaturePoints(); auto tmp_side_coords_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),tmp_side_coords.get_view()); - TEST_FLOATING_EQUALITY(tmp_side_coords_host(0,0,0), + TEST_FLOATING_EQUALITY(tmp_side_coords_host(0,0,0), realspace_x_coord_2, 1.0e-8); - TEST_FLOATING_EQUALITY(tmp_side_coords_host(0,0,1), + TEST_FLOATING_EQUALITY(tmp_side_coords_host(0,0,1), realspace_y_coord_2, 1.0e-8); } TEUCHOS_UNIT_TEST(integration_values, control_volume_boundary) { - Teuchos::RCP topo = + Teuchos::RCP topo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >())); const int num_cells = 2; @@ -211,16 +205,13 @@ namespace panzer { const panzer::CellData cell_data(num_cells,cell_side,topo); std::string cv_type = "boundary"; - RCP int_rule_bc = + RCP int_rule_bc = rcp(new IntegrationRule(cell_data, cv_type)); - panzer::IntegrationValues2 int_values_bc("prefix_",true); panzer::MDFieldArrayFactory af("prefix_",true); - int_values_bc.setupArrays(int_rule_bc); - const int num_vertices = int_rule_bc->topology->getNodeCount(); - PHX::MDField node_coordinates + PHX::MDField node_coordinates = af.buildStaticArray("nc",num_cells, num_vertices, base_cell_dimension); // Set up node coordinates. Here we assume the following @@ -248,187 +239,32 @@ namespace panzer { node_coordinates_k(cell,3,y) = 1.0; }); - int_values_bc.evaluateValues(node_coordinates); - auto ip_coords_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),int_values_bc.ip_coordinates.get_view()); - - TEST_EQUALITY(int_values_bc.ip_coordinates.extent(1), 2); + panzer::IntegrationValues2 int_values_bc("prefix_"); + int_values_bc.setup(int_rule_bc, node_coordinates); + + auto ip_coords_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),int_values_bc.getCubaturePoints().get_view()); + + out << "Points:\n"; + for(int i=0; i topo = - Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >())); - - const int in_num_cells = 20; - const int base_cell_dimension = 2; - const panzer::CellData cell_data(in_num_cells,topo); - - const int cubature_degree = 2; - RCP int_rule = - rcp(new IntegrationRule(cubature_degree, cell_data)); - - panzer::IntegrationValues2 int_values("prefix_",true); - panzer::MDFieldArrayFactory af("prefix_",true); - - int_values.setupArrays(int_rule); - - const int num_vertices = int_rule->topology->getNodeCount(); - PHX::MDField node_coordinates - = af.buildStaticArray("nc",in_num_cells, num_vertices, base_cell_dimension); - - // Set up node coordinates. Here we assume the following - // ordering. This needs to be consistent with shards topology, - // otherwise we will get negative determinates - - // 3(0,1)---2(1.5,1.5) - // | 0 | - // | | - // 0(0,0)---1(1,0) - - typedef panzer::ArrayTraits >::size_type size_type; - const size_type x = 0; - const size_type y = 1; - auto node_coordinates_k = node_coordinates.get_view(); - Kokkos::parallel_for("initialize node coords",node_coordinates.extent(0), - KOKKOS_LAMBDA (const int cell) { - node_coordinates_k(cell,0,x) = 0.0; - node_coordinates_k(cell,0,y) = 0.0; - node_coordinates_k(cell,1,x) = 1.0; - node_coordinates_k(cell,1,y) = 0.0; - node_coordinates_k(cell,2,x) = 1.5; - node_coordinates_k(cell,2,y) = 1.5; - node_coordinates_k(cell,3,x) = 0.0; - node_coordinates_k(cell,3,y) = 1.0; - }); - - int_values.evaluateValues(node_coordinates); - - // pull out arrays to look at - - auto weighted_measure = Kokkos::create_mirror_view(int_values.weighted_measure.get_view()); - auto jac_det = Kokkos::create_mirror_view(int_values.jac_det.get_view()); - auto ref_ip_coordinates = Kokkos::create_mirror_view(int_values.ref_ip_coordinates.get_view()); - auto ip_coordinates = Kokkos::create_mirror_view(int_values.ip_coordinates.get_view()); - auto jac = Kokkos::create_mirror_view(int_values.jac.get_view()); - auto jac_inv = Kokkos::create_mirror_view(int_values.jac_inv.get_view()); - - int num_cells = ip_coordinates.extent(0); - int num_points = ip_coordinates.extent(1); - int num_dim = ip_coordinates.extent(2); - - TEST_EQUALITY(num_cells,in_num_cells); - TEST_EQUALITY(num_points,4); - TEST_EQUALITY(num_dim,2); - - int ref_cell = 0; - - { - int tst_cell = 1; - int org_pt = 0; - int new_pt = 1; - - auto ref_ip_coordinates_k = int_values.ref_ip_coordinates.get_view(); - auto ip_coordinates_k = int_values.ip_coordinates.get_view(); - auto weighted_measure_k = int_values.weighted_measure.get_view(); - auto jac_k = int_values.jac.get_view(); - auto jac_det_k = int_values.jac_det.get_view(); - auto jac_inv_k = int_values.jac_inv.get_view(); - auto surface_normals_k = int_values.surface_normals.get_view(); - auto surface_rotation_matrices_k = int_values.surface_rotation_matrices.get_view(); - Kokkos::parallel_for("swap qp",1,KOKKOS_LAMBDA (const int ) { - panzer::swapQuadraturePoints(tst_cell,org_pt,new_pt, - ref_ip_coordinates_k, - ip_coordinates_k, - weighted_measure_k, - jac_k, - jac_det_k, - jac_inv_k, - surface_normals_k, - surface_rotation_matrices_k); - }); - - Kokkos::deep_copy(weighted_measure,int_values.weighted_measure.get_view()); - Kokkos::deep_copy(jac_det,int_values.jac_det.get_view()); - Kokkos::deep_copy(ref_ip_coordinates,int_values.ref_ip_coordinates.get_view()); - Kokkos::deep_copy(ip_coordinates,int_values.ip_coordinates.get_view()); - Kokkos::deep_copy(jac,int_values.jac.get_view()); - Kokkos::deep_copy(jac_inv,int_values.jac_inv.get_view()); - - TEST_EQUALITY( weighted_measure(ref_cell,org_pt), weighted_measure(tst_cell,new_pt)); - TEST_EQUALITY( jac_det(ref_cell,org_pt), jac_det(tst_cell,new_pt)); - TEST_EQUALITY(ref_ip_coordinates(ref_cell,org_pt,0), ref_ip_coordinates(tst_cell,new_pt,0)); - TEST_EQUALITY(ref_ip_coordinates(ref_cell,org_pt,1), ref_ip_coordinates(tst_cell,new_pt,1)); - TEST_EQUALITY( ip_coordinates(ref_cell,org_pt,0), ip_coordinates(tst_cell,new_pt,0)); - TEST_EQUALITY( ip_coordinates(ref_cell,org_pt,1), ip_coordinates(tst_cell,new_pt,1)); - TEST_EQUALITY( jac(ref_cell,org_pt,0,0), jac(tst_cell,new_pt,0,0)); - TEST_EQUALITY( jac(ref_cell,org_pt,0,1), jac(tst_cell,new_pt,0,1)); - TEST_EQUALITY( jac(ref_cell,org_pt,1,0), jac(tst_cell,new_pt,1,0)); - TEST_EQUALITY( jac(ref_cell,org_pt,1,1), jac(tst_cell,new_pt,1,1)); - TEST_EQUALITY( jac_inv(ref_cell,org_pt,0,0), jac_inv(tst_cell,new_pt,0,0)); - TEST_EQUALITY( jac_inv(ref_cell,org_pt,0,1), jac_inv(tst_cell,new_pt,0,1)); - TEST_EQUALITY( jac_inv(ref_cell,org_pt,1,0), jac_inv(tst_cell,new_pt,1,0)); - TEST_EQUALITY( jac_inv(ref_cell,org_pt,1,1), jac_inv(tst_cell,new_pt,1,1)); - } - - { - int tst_cell = 7; - int org_pt = 1; - int new_pt = 3; - - auto ref_ip_coordinates_k = int_values.ref_ip_coordinates.get_view(); - auto ip_coordinates_k = int_values.ip_coordinates.get_view(); - auto weighted_measure_k = int_values.weighted_measure.get_view(); - auto jac_k = int_values.jac.get_view(); - auto jac_det_k = int_values.jac_det.get_view(); - auto jac_inv_k = int_values.jac_inv.get_view(); - auto surface_normals_k = int_values.surface_normals.get_view(); - auto surface_rotation_matrices_k = int_values.surface_rotation_matrices.get_view(); - // Kokkos::parallel_for("swap qp",1,DoSwapOnDevice(tst_cell,org_pt,new_pt,int_values)); - Kokkos::parallel_for("swap qp",1,KOKKOS_LAMBDA (const int ) { - panzer::swapQuadraturePoints(tst_cell,org_pt,new_pt, - ref_ip_coordinates_k, - ip_coordinates_k, - weighted_measure_k, - jac_k, - jac_det_k, - jac_inv_k, - surface_normals_k, - surface_rotation_matrices_k); - }); - - Kokkos::deep_copy(weighted_measure,int_values.weighted_measure.get_view()); - Kokkos::deep_copy(jac_det,int_values.jac_det.get_view()); - Kokkos::deep_copy(ref_ip_coordinates,int_values.ref_ip_coordinates.get_view()); - Kokkos::deep_copy(ip_coordinates,int_values.ip_coordinates.get_view()); - Kokkos::deep_copy(jac,int_values.jac.get_view()); - Kokkos::deep_copy(jac_inv,int_values.jac_inv.get_view()); - - - TEST_EQUALITY( weighted_measure(ref_cell,org_pt), weighted_measure(tst_cell,new_pt)); - TEST_EQUALITY( jac_det(ref_cell,org_pt), jac_det(tst_cell,new_pt)); - TEST_EQUALITY(ref_ip_coordinates(ref_cell,org_pt,0), ref_ip_coordinates(tst_cell,new_pt,0)); - TEST_EQUALITY(ref_ip_coordinates(ref_cell,org_pt,1), ref_ip_coordinates(tst_cell,new_pt,1)); - TEST_EQUALITY( ip_coordinates(ref_cell,org_pt,0), ip_coordinates(tst_cell,new_pt,0)); - TEST_EQUALITY( ip_coordinates(ref_cell,org_pt,1), ip_coordinates(tst_cell,new_pt,1)); - TEST_EQUALITY( jac(ref_cell,org_pt,0,0), jac(tst_cell,new_pt,0,0)); - TEST_EQUALITY( jac(ref_cell,org_pt,0,1), jac(tst_cell,new_pt,0,1)); - TEST_EQUALITY( jac(ref_cell,org_pt,1,0), jac(tst_cell,new_pt,1,0)); - TEST_EQUALITY( jac(ref_cell,org_pt,1,1), jac(tst_cell,new_pt,1,1)); - TEST_EQUALITY( jac_inv(ref_cell,org_pt,0,0), jac_inv(tst_cell,new_pt,0,0)); - TEST_EQUALITY( jac_inv(ref_cell,org_pt,0,1), jac_inv(tst_cell,new_pt,0,1)); - TEST_EQUALITY( jac_inv(ref_cell,org_pt,1,0), jac_inv(tst_cell,new_pt,1,0)); - TEST_EQUALITY( jac_inv(ref_cell,org_pt,1,1), jac_inv(tst_cell,new_pt,1,1)); - } - } - TEUCHOS_UNIT_TEST(integration_values, surface_quadrature){ /* @@ -440,7 +276,7 @@ namespace panzer { \ | / \ | / 3--4--5 - + Cell 0: 0,3,4,1 Cell 1: 1,4,5,2 @@ -452,7 +288,7 @@ namespace panzer { 3: 0, -1 (1,0) 4: 1, -1 (4,5) 5: 1, -1 (2,1) - + */ // First we build the mesh @@ -518,11 +354,9 @@ namespace panzer { cell_to_faces_host(1,0) = 2; cell_to_faces_host(1,1) = 4; cell_to_faces_host(1,2) = 0; cell_to_faces_host(1,3) = 5; Kokkos::deep_copy(mesh.cell_to_faces,cell_to_faces_host); } - + const auto id = panzer::IntegrationDescriptor(2, panzer::IntegrationDescriptor::SURFACE); auto int_rule = Teuchos::rcp(new IntegrationRule(id, mesh.cell_topology, 2, 6)); - panzer::IntegrationValues2 int_values("prefix_",true); - int_values.setupArrays(int_rule); // Fill node coordinates panzer::MDFieldArrayFactory af("prefix_",true); @@ -537,6 +371,9 @@ namespace panzer { Kokkos::deep_copy(node_coordinates.get_static_view(),node_coordinates_host); } + panzer::IntegrationValues2 int_values("prefix_"); + int_values.setup(int_rule, node_coordinates); + // Make sure the integration values fail to construct without the subcell connectivity TEST_THROW(int_values.evaluateValues(node_coordinates), std::logic_error); @@ -545,11 +382,11 @@ namespace panzer { connectivity->setup(mesh); // Now we try again - int_values.evaluateValues(node_coordinates,-1,connectivity); - + int_values.evaluateValues(node_coordinates,-1,connectivity,0); + // This test will focus on normals and points - const auto normals = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),int_values.surface_normals.get_static_view()); - const auto points = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),int_values.ip_coordinates.get_static_view()); + const auto normals = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),int_values.getSurfaceNormals().get_static_view()); + const auto points = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),int_values.getCubaturePoints().get_static_view()); const int num_points_per_cell = points.extent_int(1); const int num_points_per_face = num_points_per_cell / 4; @@ -568,7 +405,7 @@ namespace panzer { // We need to make sure that face 0 and 2 have aligned quadrature points const double tolerance = 1.e-14; - + // Face 0 { const int cell_0 = connectivity->cellForSubcellHost(0,0); @@ -588,23 +425,26 @@ namespace panzer { const double normal_1[2] = { 2./sqrt5,-1./sqrt5}; // Note that y will be equal, but x will be different - + out << "Comparing cell_0 "< 0.5); TEST_FLOATING_EQUALITY(points(cell_0,point_0,1), points(cell_1,point_1,1), tolerance); - TEST_FLOATING_EQUALITY(normals(cell_0,point_0,0), normal_0[0], tolerance); + TEST_FLOATING_EQUALITY(normals(cell_0,point_0,0), normal_0[0], tolerance); TEST_FLOATING_EQUALITY(normals(cell_0,point_0,1), normal_0[1], tolerance); - TEST_FLOATING_EQUALITY(normals(cell_1,point_1,0), normal_1[0], tolerance); + TEST_FLOATING_EQUALITY(normals(cell_1,point_1,0), normal_1[0], tolerance); TEST_FLOATING_EQUALITY(normals(cell_1,point_1,1), normal_1[1], tolerance); } } - + // Face 2 { const int cell_0 = connectivity->cellForSubcellHost(2,0); @@ -621,17 +461,21 @@ namespace panzer { const double normal_0[2] = { 1.,0.}; const double normal_1[2] = {-1.,0.}; - + + out << "Comparing cell_0 "<