From 0f8e903fd37b9f7a3f01f49b31234d48ef69f16e Mon Sep 17 00:00:00 2001 From: "Eric C. Cyr" Date: Fri, 15 Jun 2018 12:05:29 -0600 Subject: [PATCH] Panzer: Fixing integration values side qp construction These were not correctly being reordered. This change fixes this. 1) Adding reorder function to utilities This function automates the reordering of an array by a specified vector 2) side points in IV2, handling ref_ip_coords in a consistent way 3) Extracting sorting function in IntegrationValues Adding unit test for 1d sorting 4) Core, fixing globbing issue Toucing CMakeLists.txt 5) Adding 2d coordinates test to IV2 This is a copy of the 1D one (but in 2D). It excercises the periodic case 6) Updating unit tests associated with sorting face quadrature rules 7) Adding swap function for quadrature points Incorporated into the reordering routines and added unit tests --- .../panzer_workset_builder/CMakeLists.txt | 7 + .../check_sorted_face_quads.cpp | 246 ++++++++++++++++++ packages/panzer/core/src/CMakeLists.txt | 1 + .../panzer/core/src/Panzer_UtilityAlgs.cpp | 30 +++ .../panzer/core/src/Panzer_UtilityAlgs.hpp | 22 ++ packages/panzer/core/test/CMakeLists.txt | 6 + packages/panzer/core/test/utility_algs.cpp | 94 +++++++ .../disc-fe/src/Panzer_IntegrationValues2.cpp | 121 +++++---- .../disc-fe/src/Panzer_IntegrationValues2.hpp | 27 ++ .../test/core_tests/integration_values2.cpp | 230 ++++++++++++++++ 10 files changed, 735 insertions(+), 49 deletions(-) create mode 100644 packages/panzer/adapters-stk/test/panzer_workset_builder/check_sorted_face_quads.cpp create mode 100644 packages/panzer/core/src/Panzer_UtilityAlgs.cpp create mode 100644 packages/panzer/core/src/Panzer_UtilityAlgs.hpp create mode 100644 packages/panzer/core/test/utility_algs.cpp diff --git a/packages/panzer/adapters-stk/test/panzer_workset_builder/CMakeLists.txt b/packages/panzer/adapters-stk/test/panzer_workset_builder/CMakeLists.txt index ba1ae8df0092..962806673338 100644 --- a/packages/panzer/adapters-stk/test/panzer_workset_builder/CMakeLists.txt +++ b/packages/panzer/adapters-stk/test/panzer_workset_builder/CMakeLists.txt @@ -54,3 +54,10 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST( COMM serial mpi NUM_MPI_PROCS 1 ) + +TRIBITS_ADD_EXECUTABLE_AND_TEST( + check_sorted_face_quads + SOURCES check_sorted_face_quads.cpp ${UNIT_TEST_DRIVER} + COMM serial mpi + NUM_MPI_PROCS 1 + ) diff --git a/packages/panzer/adapters-stk/test/panzer_workset_builder/check_sorted_face_quads.cpp b/packages/panzer/adapters-stk/test/panzer_workset_builder/check_sorted_face_quads.cpp new file mode 100644 index 000000000000..452c2a258456 --- /dev/null +++ b/packages/panzer/adapters-stk/test/panzer_workset_builder/check_sorted_face_quads.cpp @@ -0,0 +1,246 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER +// + +#include + +#include +#include +#include +#include + +using Teuchos::RCP; +using Teuchos::rcp; + +#include "Teuchos_DefaultComm.hpp" +#include "Teuchos_GlobalMPISession.hpp" + +#include "Panzer_STK_Version.hpp" +#include "PanzerAdaptersSTK_config.hpp" + +#include "Panzer_WorksetContainer.hpp" +#include "Panzer_WorksetDescriptor.hpp" +#include "Panzer_WorksetNeeds.hpp" +#include "Panzer_IntrepidBasisFactory.hpp" +#include "Panzer_IntegrationDescriptor.hpp" +#include "Panzer_DOFManager.hpp" +#include "Panzer_SubcellConnectivity.hpp" + + +#include "Panzer_STK_Interface.hpp" +#include "Panzer_STK_SquareQuadMeshFactory.hpp" +#include "Panzer_STKConnManager.hpp" +#include "Panzer_STK_WorksetFactory.hpp" + +namespace panzer { + + TEUCHOS_UNIT_TEST(check_sorted_face_quads, basic) + { + using Teuchos::RCP; + using Teuchos::rcp; + + typedef Intrepid2::Basis IntrepidBasis; + + std::string element_block = "eblock-0_0"; + std::string sideset = "left"; + + RCP pl = rcp(new Teuchos::ParameterList); + pl->set ("X Blocks", 1); + pl->set ("X Elements", 128); + pl->set("X0", 0.0); + pl->set("Xf", 2.5); + pl->set ("X Procs", 1); + pl->set ("Y Blocks", 1); + pl->set ("Y Elements", 2); + pl->set("Y0", 0); + pl->set("Yf", 1.11804); + pl->set ("Y Procs", 1); + pl->sublist("Periodic BCs"); + pl->sublist("Periodic BCs").set("Periodic Condition 1","y-all 1e-8: left;right"); + pl->sublist("Periodic BCs").set("Periodic Condition 2","x-all 1e-8: top;bottom"); + pl->sublist("Periodic BCs").set("Count",2); + + panzer_stk::SquareQuadMeshFactory factory; + factory.setParameterList(pl); + RCP mesh = factory.buildMesh(MPI_COMM_WORLD); + RCP > tComm = Teuchos::rcp(new Teuchos::MpiComm(MPI_COMM_WORLD)); + + // build DOF Manager (with a single HDiv basis) + ///////////////////////////////////////////////////////////// + + // build the connection manager + const RCP > + conn_manager = rcp(new panzer_stk::STKConnManager(mesh)); + + RCP > dof_manager + = rcp(new panzer::DOFManager(conn_manager,MPI_COMM_WORLD)); + + // build an intrepid basis and a related field pattern for seeding the DOFManager + { + RCP intrepid_basis + = panzer::createIntrepid2Basis("HGrad",1, + *mesh->getCellTopology(element_block)); + RCP field_pattern = rcp(new panzer::Intrepid2FieldPattern(intrepid_basis)); + + dof_manager->addField(element_block, "T", field_pattern); + } + dof_manager->buildGlobalUnknowns(); + + // build WorksetContainer + ////////////////////////////////////////////////////////////// + + panzer::IntegrationDescriptor sid(3, panzer::IntegrationDescriptor::SURFACE); + panzer::BasisDescriptor bd(2, "HGrad"); + std::map wkstRequirements; + wkstRequirements[element_block].addIntegrator(sid); + wkstRequirements[element_block].addBasis(bd); + + 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)); + + wkstContainer->setGlobalIndexer(dof_manager); + + panzer::WorksetDescriptor workset_descriptor(element_block, panzer::WorksetSizeType::ALL_ELEMENTS, true,false); + + auto worksets = wkstContainer->getWorksets(workset_descriptor); + + TEST_ASSERT(worksets->size()==1); + + auto & workset = (*worksets)[0]; + auto rot_matrices = workset.getIntegrationValues(sid).surface_rotation_matrices; + auto normals = workset.getIntegrationValues(sid).surface_normals; + auto ip_coordinates = workset.getIntegrationValues(sid).ip_coordinates; + + const panzer::SubcellConnectivity & face_connectivity = workset.getFaceConnectivity(); + + + const int num_points = normals.extent(1); + const int num_faces = face_connectivity.numSubcells(); + const int num_faces_per_cell = face_connectivity.numSubcellsOnCell(0); + const int num_points_per_face = num_points / num_faces_per_cell; + + TEST_EQUALITY(num_faces,512); // 128 * 2 vertical edges (2 removed by periodicity) + // 128 * 2 horizontal edges (128 removed by periodicty) + TEST_EQUALITY(num_faces_per_cell,4); + TEST_EQUALITY(num_points,2*4); + + for(int face=0; face & order,std::function swapper) +{ + // each entry has to be sorted + for(std::size_t i=0;i +#include + +namespace panzer{ + +/** + * \brief Using a functor, reorder an array using a order vector. + * + * \param[in,out] order Vector that describes the desired order. + * Note on output, this array will be sorted. + * \param[in] swapper Functor that swaps entries in the array to be + * sorted. Data being swapped must be captured + * by the functor. + */ +void reorder(std::vector & order,std::function swapper); + +} + +#endif diff --git a/packages/panzer/core/test/CMakeLists.txt b/packages/panzer/core/test/CMakeLists.txt index 3a7cfb4022ff..b7cea56901db 100644 --- a/packages/panzer/core/test/CMakeLists.txt +++ b/packages/panzer/core/test/CMakeLists.txt @@ -24,3 +24,9 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST( NUM_MPI_PROCS 1 ) +TRIBITS_ADD_EXECUTABLE_AND_TEST( + utility_algs + SOURCES utility_algs.cpp ${UNIT_TEST_DRIVER} + NUM_MPI_PROCS 1 + ) + diff --git a/packages/panzer/core/test/utility_algs.cpp b/packages/panzer/core/test/utility_algs.cpp new file mode 100644 index 000000000000..837b686f33ed --- /dev/null +++ b/packages/panzer/core/test/utility_algs.cpp @@ -0,0 +1,94 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#include +#include +#include "Panzer_UtilityAlgs.hpp" + +namespace panzer { + + TEUCHOS_UNIT_TEST(utility_algs, reorder) + { + // test 1 + { + std::vector order = {2, 0, 1}; + std::vector data = {'a', 'b', 'c'}; + + // run reorder algorithm + reorder(order,[&data](int a,int b) { + std::swap(data[a],data[b]); + }); + + TEST_ASSERT(order==std::vector({0,1,2})); + TEST_ASSERT(data==std::vector({'c','a','b'})); + } + + // test 2 + { + std::vector order = {0, 1}; + std::vector data = {'b', 'c'}; + + // run reorder algorithm + reorder(order,[&data](int a,int b) { + std::swap(data[a],data[b]); + }); + + TEST_ASSERT(order==std::vector({0,1})); + TEST_ASSERT(data==std::vector({'b','c'})); + } + + // test 3 + { + std::vector order = {1, 0}; + std::vector data = {'b', 'c'}; + + // run reorder algorithm + reorder(order,[&data](int a,int b) { + std::swap(data[a],data[b]); + }); + + TEST_ASSERT(order==std::vector({0,1})); + TEST_ASSERT(data==std::vector({'c','b'})); + } + } + +} // end namespace panzer diff --git a/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.cpp b/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.cpp index f30302e6971e..fe46aff6423a 100644 --- a/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.cpp +++ b/packages/panzer/disc-fe/src/Panzer_IntegrationValues2.cpp @@ -41,6 +41,7 @@ // @HEADER #include "Panzer_IntegrationValues2.hpp" +#include "Panzer_UtilityAlgs.hpp" #include "Shards_CellTopology.hpp" @@ -495,6 +496,71 @@ convertNormalToRotationMatrix(const T normal[3], T transverse[3], T binormal[3]) } +template +void IntegrationValues2:: +swapQuadraturePoints(int cell, + int a, + int b) +{ + const int new_cell_point = a; + const int old_cell_point = b; + + const int cell_dim = ref_ip_coordinates.extent(2); + + 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 +void IntegrationValues2:: +uniqueCoordOrdering(Array_CellIPDim & coords, + int cell, + int offset, + std::vector & order) +{ + for(size_t point_index=0;point_index sorter(coords,cell,offset); + std::sort(order.begin(),order.end(),sorter); +} template void IntegrationValues2:: @@ -660,8 +726,9 @@ generateSurfaceCubatureValues(const PHX::MDField& in_node_ weighted_measure(cell,cell_point) = side_weighted_measure(cell,side_point); jac_det(cell,cell_point) = side_det_jacobian(cell,side_point); for(int dim=0;dim& in_node_ std::vector point_indexes(num_points_on_face,-1); for(int cell=0; cell sorter(ip_coordinates,cell,point_offset); - std::sort(point_indexes.begin(),point_indexes.end(),sorter); + // build a point index array based on point coordinates + uniqueCoordOrdering(ip_coordinates,cell,point_offset,point_indexes); // Indexes are now sorted, now we swap everything around - for(int old_point_idx=0;old_point_idx size: span() == jac.span() + + /** + * \brief Using coordinate build an arrray that specifies a unique ordering. + * + * Used for side integration points. Compute a unique ordering in a cell and + * point offset. + * + * \param[in] coords Coordinates array (cell,IP,Dim) + * \param[in] cell Cell index + * \param[in] offset Offset into the points + * \param[out] order Ordering array on output, correctly sized on input + * (offset + order.size() <= coords.extent(1)) + */ + static void uniqueCoordOrdering(Array_CellIPDim & coords, + int cell, + int offset, + std::vector & order); + + /** + * \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 + */ + void swapQuadraturePoints(int cell,int a,int b); + protected: 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 1ba8162b1d61..ff53747afebf 100644 --- a/packages/panzer/disc-fe/test/core_tests/integration_values2.cpp +++ b/packages/panzer/disc-fe/test/core_tests/integration_values2.cpp @@ -245,4 +245,234 @@ namespace panzer { realspace_y_coord_1, 1.0e-8); } + + TEUCHOS_UNIT_TEST(integration_values, coord_ordering_1d) + { + typedef IntegrationValues2 IV; + MDFieldArrayFactory af("",true); + + int num_faces = 2; + int num_points_per_face = 4; + std::vector order(num_points_per_face,-1); + + IV::Array_CellIPDim coords = af.template buildStaticArray("coord",2,num_faces*num_points_per_face,1); + coords(0,0,0) = -1.0; coords(0,1,0) = 0.0; coords(0,2,0) = 1.0; coords(0,3,0) = 2.0; + coords(0,4,0) = 1.0; coords(0,5,0) = 2.0; coords(0,6,0) = -1.0; coords(0,7,0) = 0.0; + coords(1,0,0) = -1.0; coords(1,1,0) = 0.0; coords(1,2,0) = 1.0; coords(1,3,0) = 2.0; + coords(1,4,0) = 1.0; coords(1,5,0) = 2.0; coords(1,6,0) = -1.0; coords(1,7,0) = 0.0; + + int cell = 0, offset = 0; + + cell = 0; offset = 0; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({0,1,2,3})); + + cell = 0; offset = 4; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({2,3,0,1})); + + cell = 1; offset = 0; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({0,1,2,3})); + + cell = 1; offset = 4; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({2,3,0,1})); + } + + TEUCHOS_UNIT_TEST(integration_values, coord_ordering_2d) + { + typedef IntegrationValues2 IV; + MDFieldArrayFactory af("",true); + + { + int num_faces = 2; + int num_points_per_face = 4; + std::vector order(num_points_per_face,-1); + + IV::Array_CellIPDim coords = af.template buildStaticArray("coord",2,num_faces*num_points_per_face,2); + + // cell 0 + coords(0,0,0) = -1.0; coords(0,1,0) = 0.0; coords(0,2,0) = 1.0; coords(0,3,0) = 2.0; + coords(0,0,1) = 0.0; coords(0,1,1) = 0.0; coords(0,2,1) = 0.0; coords(0,3,1) = 0.0; + + coords(0,4,0) = 1.0; coords(0,5,0) = 2.0; coords(0,6,0) = -1.0; coords(0,7,0) = 0.0; + coords(0,4,1) = 0.0; coords(0,5,1) = 0.0; coords(0,6,1) = 0.0; coords(0,7,1) = 0.0; + + // cell 1 + coords(1,0,0) = 2.0; coords(1,1,0) = 2.0; coords(1,2,0) = 2.0; coords(1,3,0) = 2.0; + coords(1,0,0) = -1.1; coords(1,1,1) = 0.0; coords(1,2,1) = 1.0; coords(1,3,1) = 2.0; + + coords(1,4,0) = 2.0; coords(1,5,0) = 2.0; coords(1,6,0) = 2.0; coords(1,7,0) = 2.0; + coords(1,4,1) = 1.0; coords(1,5,1) = 2.0; coords(1,6,1) = -1.0; coords(1,7,1) = 0.0; + + int cell = 0, offset = 0; + + cell = 0; offset = 0; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({0,1,2,3})); + + cell = 0; offset = 4; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({2,3,0,1})); + + cell = 1; offset = 0; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({0,1,2,3})); + + cell = 1; offset = 4; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({2,3,0,1})); + } + + // this was the original failing case + { + std::vector order(2,-1); + + IV::Array_CellIPDim coords = af.template buildStaticArray("coord",4,2,2); + + coords(0,0,0) = 0.0; coords(0,0,1) = 4.4088517374119224e-01; + coords(0,1,0) = 0.0; coords(0,1,1) = 1.1813482625880771e-01; + + coords(1,0,0) = 2.5; coords(1,0,1) = 4.4088517374119224e-01; + coords(1,1,0) = 2.5; coords(1,1,1) = 1.1813482625880771e-01; + + coords(2,0,0) = 0.0; coords(2,0,1) = 1.1813482625880771e-01; + coords(2,1,0) = 0.0; coords(2,1,1) = 4.4088517374119224e-01; + + coords(3,0,0) = 2.5; coords(3,0,1) = 1.1813482625880771e-01; + coords(3,1,0) = 2.5; coords(3,1,1) = 4.4088517374119224e-01; + + int cell = 0, offset = 0; + + cell = 0; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({1,0})); + + cell = 1; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({1,0})); + + cell = 2; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({0,1})); + + cell = 3; + IV::uniqueCoordOrdering(coords,cell,offset,order); + TEST_ASSERT(order==std::vector({0,1})); + } + } + + TEUCHOS_UNIT_TEST(integration_values, quadpt_swap) + { + Teuchos::RCP 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; + for (size_type cell = 0; cell < node_coordinates.extent(0); ++cell) { + node_coordinates(cell,0,x) = 0.0; + node_coordinates(cell,0,y) = 0.0; + node_coordinates(cell,1,x) = 1.0; + node_coordinates(cell,1,y) = 0.0; + node_coordinates(cell,2,x) = 1.5; + node_coordinates(cell,2,y) = 1.5; + node_coordinates(cell,3,x) = 0.0; + node_coordinates(cell,3,y) = 1.0; + } + + int_values.evaluateValues(node_coordinates); + + // pull out arrays to look at + + auto weighted_measure = int_values.weighted_measure; + auto jac_det = int_values.jac_det; + auto ref_ip_coordinates = int_values.ref_ip_coordinates; + auto ip_coordinates = int_values.ip_coordinates; + auto jac = int_values.jac; + auto jac_inv = int_values.jac_inv; + + 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; + + int_values.swapQuadraturePoints(tst_cell,org_pt,new_pt); + + 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; + + int_values.swapQuadraturePoints(tst_cell,org_pt,new_pt); + + 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)); + } + } }