From 13092ecb26fb469ac088b42b9553e50face1f4e3 Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Thu, 5 Sep 2024 14:06:45 -0600 Subject: [PATCH] Intrepid2: Fix inclusion test Fixed some issues exposed by a Sierra test using float instead of double. Added Cell inclusion tests using the float type. --- .../src/Cell/Intrepid2_CellDataDef.hpp | 31 ++- .../Cell/Intrepid2_CellToolsDefInclusion.hpp | 12 +- ...id2_CubatureDirectLineGaussJacobi20Def.hpp | 2 +- .../Intrepid2_CubatureDirectTetDefaultDef.hpp | 2 +- ...ntrepid2_CubatureDirectTetSymmetricDef.hpp | 2 +- .../Intrepid2_CubatureDirectTriDefaultDef.hpp | 2 +- ...ntrepid2_CubatureDirectTriSymmetricDef.hpp | 2 +- .../intrepid2/src/Shared/Intrepid2_Types.hpp | 2 +- .../intrepid2/unit-test/Cell/CMakeLists.txt | 4 + packages/intrepid2/unit-test/Cell/test_01.hpp | 2 +- packages/intrepid2/unit-test/Cell/test_02.hpp | 4 +- packages/intrepid2/unit-test/Cell/test_03.hpp | 2 +- packages/intrepid2/unit-test/Cell/test_04.hpp | 4 +- packages/intrepid2/unit-test/Cell/test_06.hpp | 2 +- packages/intrepid2/unit-test/Cell/test_07.hpp | 247 +++++++++--------- 15 files changed, 164 insertions(+), 156 deletions(-) diff --git a/packages/intrepid2/src/Cell/Intrepid2_CellDataDef.hpp b/packages/intrepid2/src/Cell/Intrepid2_CellDataDef.hpp index d21e77553fac..7e14c0e9bb4b 100644 --- a/packages/intrepid2/src/Cell/Intrepid2_CellDataDef.hpp +++ b/packages/intrepid2/src/Cell/Intrepid2_CellDataDef.hpp @@ -831,6 +831,7 @@ refCenterDataStatic_ = { bool PointInclusion::key>:: check(const PointViewType &point, const ScalarType threshold) { + //this implementation should work when PointType is a Sacado Fad const ScalarType minus_one = -1.0 - threshold, plus_one = 1.0 + threshold; return (minus_one <= point(0) && point(0) <= plus_one); } @@ -840,7 +841,10 @@ refCenterDataStatic_ = { bool PointInclusion::key>:: check(const PointViewType &point, const ScalarType threshold) { - const ScalarType distance = max( max( -point(0), -point(1) ), ScalarType( point(0) + point(1) - 1.0 ) ); + //this implementation should work when PointType is a Sacado Fad + using PointType = typename PointViewType::value_type; + const PointType one = 1.0; + const PointType distance = max( max( -point(0), -point(1) ), point(0) + point(1) - one ); return distance < threshold; } @@ -848,8 +852,8 @@ refCenterDataStatic_ = { KOKKOS_INLINE_FUNCTION bool PointInclusion::key>:: - check(const PointViewType &point, - const ScalarType threshold) { + check(const PointViewType &point, const ScalarType threshold) { + //this implementation should work when PointType is a Sacado Fad const ScalarType minus_one = -1.0 - threshold, plus_one = 1.0 + threshold; return ((minus_one <= point(0) && point(0) <= plus_one) && (minus_one <= point(1) && point(1) <= plus_one)); @@ -860,8 +864,11 @@ refCenterDataStatic_ = { bool PointInclusion::key>:: check(const PointViewType &point, const ScalarType threshold) { - const ScalarType distance = max( max(-point(0),-point(1)), - max(-point(2), point(0) + point(1) + point(2) - 1) ); + //this implementation should work when PointType is a Sacado Fad + using PointType = typename PointViewType::value_type; + const PointType one = 1.0; + const PointType distance = max( max(-point(0),-point(1)), + max(-point(2), point(0) + point(1) + point(2) - one) ); return distance < threshold; } @@ -870,6 +877,7 @@ refCenterDataStatic_ = { bool PointInclusion::key>:: check(const PointViewType &point, const ScalarType threshold) { + //this implementation should work when PointType is a Sacado Fad const ScalarType minus_one = -1.0 - threshold, plus_one = 1.0 + threshold; return ((minus_one <= point(0) && point(0) <= plus_one) && (minus_one <= point(1) && point(1) <= plus_one) && @@ -881,9 +889,11 @@ refCenterDataStatic_ = { bool PointInclusion::key>:: check(const PointViewType &point, const ScalarType threshold) { - const ScalarType minus_one = -1.0 - threshold, plus_one = 1.0 + threshold, minus_zero = -threshold; - const ScalarType left = minus_one + point(2); - const ScalarType right = plus_one - point(2); + //this implementation should work when PointType is a Sacado Fad + using PointType = typename PointViewType::value_type; + const PointType minus_one = -1.0 - threshold, plus_one = 1.0 + threshold, minus_zero = -threshold; + const PointType left = minus_one + point(2); + const PointType right = plus_one - point(2); return ((left <= point(0) && point(0) <= right) && (left <= point(1) && point(1) <= right) && (minus_zero <= point(2) && point(2) <= plus_one)); @@ -894,8 +904,11 @@ refCenterDataStatic_ = { bool PointInclusion::key>:: check(const PointViewType &point, const ScalarType threshold) { + //this implementation should work when PointType is a Sacado Fad + using PointType = typename PointViewType::value_type; const ScalarType minus_one = -1.0 - threshold, plus_one = 1.0 + threshold; - const ScalarType distance = max( max( -point(0), -point(1) ), point(0) + point(1) - 1 ); + const PointType one = 1.0; + const PointType distance = max( max( -point(0), -point(1) ), point(0) + point(1) - one ); return (distance < threshold && (minus_one <= point(2) && point(2) <= plus_one)); } diff --git a/packages/intrepid2/src/Cell/Intrepid2_CellToolsDefInclusion.hpp b/packages/intrepid2/src/Cell/Intrepid2_CellToolsDefInclusion.hpp index 1d9ecfe94b63..dfb3e04b0cf0 100644 --- a/packages/intrepid2/src/Cell/Intrepid2_CellToolsDefInclusion.hpp +++ b/packages/intrepid2/src/Cell/Intrepid2_CellToolsDefInclusion.hpp @@ -184,6 +184,10 @@ namespace Intrepid2 { checkPointwiseInclusion::key,decltype(inCell),decltype(points)>(inCell, points, threshold); break; + case shards::Tetrahedron<>::key : + checkPointwiseInclusion::key,decltype(inCell),decltype(points)>(inCell, points, threshold); + break; + case shards::Hexahedron<>::key : checkPointwiseInclusion::key,decltype(inCell),decltype(points)>(inCell, points, threshold); break; @@ -197,13 +201,7 @@ namespace Intrepid2 { break; default: - INTREPID2_TEST_FOR_EXCEPTION( !( (key == shards::Line<>::key ) || - (key == shards::Triangle<>::key) || - (key == shards::Quadrilateral<>::key) || - (key == shards::Tetrahedron<>::key) || - (key == shards::Hexahedron<>::key) || - (key == shards::Wedge<>::key) || - (key == shards::Pyramid<>::key) ), + INTREPID2_TEST_FOR_EXCEPTION( false, std::invalid_argument, ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): Invalid cell topology. "); } diff --git a/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectLineGaussJacobi20Def.hpp b/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectLineGaussJacobi20Def.hpp index 54b759bf74cc..ce8ec02ed565 100644 --- a/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectLineGaussJacobi20Def.hpp +++ b/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectLineGaussJacobi20Def.hpp @@ -18,7 +18,7 @@ namespace Intrepid2 { template CubatureDirectLineGaussJacobi20:: CubatureDirectLineGaussJacobi20(const ordinal_type degree) - : CubatureDirect
(degree, 1) { + : CubatureDirect(degree, 1) { INTREPID2_TEST_FOR_EXCEPTION( degree < 0 || degree > static_cast(Parameters::MaxCubatureDegreePyr), std::out_of_range, diff --git a/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTetDefaultDef.hpp b/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTetDefaultDef.hpp index 6a87323e8991..7a046840c9f3 100644 --- a/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTetDefaultDef.hpp +++ b/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTetDefaultDef.hpp @@ -18,7 +18,7 @@ namespace Intrepid2 { template CubatureDirectTetDefault:: CubatureDirectTetDefault(const ordinal_type degree) - : CubatureDirect
(degree, 3) { + : CubatureDirect(degree, 3) { INTREPID2_TEST_FOR_EXCEPTION( degree < 0 || degree > static_cast(Parameters::MaxCubatureDegreeTet), std::out_of_range, diff --git a/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTetSymmetricDef.hpp b/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTetSymmetricDef.hpp index 13231571a545..a13db740758a 100644 --- a/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTetSymmetricDef.hpp +++ b/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTetSymmetricDef.hpp @@ -18,7 +18,7 @@ namespace Intrepid2 { template CubatureDirectTetSymmetric:: CubatureDirectTetSymmetric(const ordinal_type degree) - : CubatureDirect
(degree, 3) { + : CubatureDirect(degree, 3) { INTREPID2_TEST_FOR_EXCEPTION( degree < 0 || degree > static_cast(Parameters::MaxCubatureDegreeTet), std::out_of_range, diff --git a/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTriDefaultDef.hpp b/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTriDefaultDef.hpp index c3c7b135a38c..378c0c44c7bf 100644 --- a/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTriDefaultDef.hpp +++ b/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTriDefaultDef.hpp @@ -18,7 +18,7 @@ namespace Intrepid2 { template CubatureDirectTriDefault:: CubatureDirectTriDefault(const ordinal_type degree) - : CubatureDirect
(degree, 2) { + : CubatureDirect(degree, 2) { INTREPID2_TEST_FOR_EXCEPTION( degree < 0 || degree >= cubatureDataStaticSize, std::out_of_range, diff --git a/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTriSymmetricDef.hpp b/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTriSymmetricDef.hpp index aa11a1212d08..0a4635738709 100644 --- a/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTriSymmetricDef.hpp +++ b/packages/intrepid2/src/Discretization/Integration/Intrepid2_CubatureDirectTriSymmetricDef.hpp @@ -18,7 +18,7 @@ namespace Intrepid2 { template CubatureDirectTriSymmetric:: CubatureDirectTriSymmetric(const ordinal_type degree) - : CubatureDirect
(degree, 2) { + : CubatureDirect(degree, 2) { INTREPID2_TEST_FOR_EXCEPTION( degree < 0 || degree >= cubatureDataStaticSize, std::out_of_range, diff --git a/packages/intrepid2/src/Shared/Intrepid2_Types.hpp b/packages/intrepid2/src/Shared/Intrepid2_Types.hpp index 8afcc665653d..8491dd36c376 100644 --- a/packages/intrepid2/src/Shared/Intrepid2_Types.hpp +++ b/packages/intrepid2/src/Shared/Intrepid2_Types.hpp @@ -56,7 +56,7 @@ namespace Intrepid2 { flt_32 s; s.f32 = 1; - s.f32++; + s.i32++; return (s.i32 < 0 ? 1 - s.f32 : s.f32 - 1); } diff --git a/packages/intrepid2/unit-test/Cell/CMakeLists.txt b/packages/intrepid2/unit-test/Cell/CMakeLists.txt index 6a72dc7540ee..a9250c464617 100644 --- a/packages/intrepid2/unit-test/Cell/CMakeLists.txt +++ b/packages/intrepid2/unit-test/Cell/CMakeLists.txt @@ -8,6 +8,10 @@ LIST(APPEND Intrepid2_TEST_ETI_VALUETYPE_NAME "DOUBLE") LIST(APPEND Intrepid2_TEST_ETI_VALUETYPE "double") LIST(APPEND Intrepid2_TEST_ETI_SACADO "0") +LIST(APPEND Intrepid2_TEST_ETI_VALUETYPE_NAME "FLOAT") +LIST(APPEND Intrepid2_TEST_ETI_VALUETYPE "float") +LIST(APPEND Intrepid2_TEST_ETI_SACADO "0") + IF (HAVE_INTREPID2_SACADO) # LIST(APPEND Intrepid2_TEST_ETI_VALUETYPE_NAME "DFAD_DOUBLE") # LIST(APPEND Intrepid2_TEST_ETI_VALUETYPE "Sacado::Fad::DFad") diff --git a/packages/intrepid2/unit-test/Cell/test_01.hpp b/packages/intrepid2/unit-test/Cell/test_01.hpp index 30a77cb71a7f..cd173cd269c8 100644 --- a/packages/intrepid2/unit-test/Cell/test_01.hpp +++ b/packages/intrepid2/unit-test/Cell/test_01.hpp @@ -200,7 +200,7 @@ namespace Intrepid2 { typedef CellTools ct; typedef Kokkos::DynRankView DynRankView; - const value_type tol = tolerence()*100.0; + const value_type tol = tolerence()*100.0; int errorFlag = 0; diff --git a/packages/intrepid2/unit-test/Cell/test_02.hpp b/packages/intrepid2/unit-test/Cell/test_02.hpp index b765e1881c46..701bcf6b96e9 100644 --- a/packages/intrepid2/unit-test/Cell/test_02.hpp +++ b/packages/intrepid2/unit-test/Cell/test_02.hpp @@ -81,7 +81,7 @@ namespace Intrepid2 { using ct = CellTools; using DynRankView = Kokkos::DynRankView; - const ValueType tol = tolerence()*100.0; + const ValueType tol = tolerence()*100.0; int errorFlag = 0; @@ -112,7 +112,7 @@ namespace Intrepid2 { { // Define cubature on the edge parametrization domain: const auto testAccuracy = 6; - CubatureDirectLineGauss edgeCubature(testAccuracy); + CubatureDirectLineGauss edgeCubature(testAccuracy); const auto cubDim = edgeCubature.getDimension(); const auto numCubPoints = edgeCubature.getNumPoints(); diff --git a/packages/intrepid2/unit-test/Cell/test_03.hpp b/packages/intrepid2/unit-test/Cell/test_03.hpp index e3bec14e9011..0170eea1d22f 100644 --- a/packages/intrepid2/unit-test/Cell/test_03.hpp +++ b/packages/intrepid2/unit-test/Cell/test_03.hpp @@ -85,7 +85,7 @@ namespace Intrepid2 { using DynRankView = Kokkos::DynRankView; using DynRankViewHost = Kokkos::DynRankView; - const ValueType tol = tolerence()*100.0; + const ValueType tol = tolerence()*100.0; int errorFlag = 0; diff --git a/packages/intrepid2/unit-test/Cell/test_04.hpp b/packages/intrepid2/unit-test/Cell/test_04.hpp index dbf8c96c68c0..b7eb9c1b61cf 100644 --- a/packages/intrepid2/unit-test/Cell/test_04.hpp +++ b/packages/intrepid2/unit-test/Cell/test_04.hpp @@ -79,7 +79,7 @@ namespace Intrepid2 { using ct = CellTools; using DynRankView = Kokkos::DynRankView; - const ValueType tol = tolerence()*100.0; + const ValueType tol = tolerence()*100.0; int errorFlag = 0; @@ -132,7 +132,7 @@ namespace Intrepid2 { continue; // 1. Define a single reference point set using cubature factory with order 4 cubature - const auto cellCubature = cubFactory.create(cell, testAccuracy); + const auto cellCubature = cubFactory.create(cell, testAccuracy); const auto cubDim = cellCubature->getDimension(); const auto cubNumPoints = cellCubature->getNumPoints(); diff --git a/packages/intrepid2/unit-test/Cell/test_06.hpp b/packages/intrepid2/unit-test/Cell/test_06.hpp index b548efc8a589..8344c0a8b5a6 100644 --- a/packages/intrepid2/unit-test/Cell/test_06.hpp +++ b/packages/intrepid2/unit-test/Cell/test_06.hpp @@ -130,7 +130,7 @@ namespace Intrepid2 { << "| |\n" << "===============================================================================\n"; - const ValueType tol = tolerence()*100.0; + const ValueType tol = tolerence()*100.0; int errorFlag = 0; diff --git a/packages/intrepid2/unit-test/Cell/test_07.hpp b/packages/intrepid2/unit-test/Cell/test_07.hpp index e6b989d32d53..813797a4475a 100644 --- a/packages/intrepid2/unit-test/Cell/test_07.hpp +++ b/packages/intrepid2/unit-test/Cell/test_07.hpp @@ -14,8 +14,6 @@ */ #include "Intrepid2_config.h" -#include "Intrepid2_CellTopologyTags.hpp" - #include "Intrepid2_DefaultCubatureFactory.hpp" #include "Teuchos_oblackholestream.hpp" @@ -37,10 +35,8 @@ namespace Intrepid2 { }; -#define INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, expected, shtopo, celltag) \ +#define INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, expected, shardsTopology) \ { \ - typedef shtopo shardsTopology; \ - typedef celltag CellTopologyTag; \ \ const auto order = 3; \ const auto cellTopo = shards::CellTopology(shards::getCellTopologyData< shardsTopology >()); \ @@ -55,14 +51,14 @@ namespace Intrepid2 { cub->getCubature(pts, wts); \ \ Kokkos::RangePolicy policy(0, P); \ - typedef F_checkPointInclusion FunctorType; \ + typedef F_checkPointInclusion FunctorType; \ Kokkos::parallel_for(policy, FunctorType(offset, check, pts)); \ \ auto check_host = Kokkos::create_mirror_view(check); \ Kokkos::deep_copy(check_host, check); \ \ for (ordinal_type i=0;i tol) { \ *outStream << "Error : checkPointInclusion at (" \ << i \ @@ -72,18 +68,19 @@ namespace Intrepid2 { } \ } \ - template + typename InputViewType> struct F_checkPointInclusion { - double _offset; + using ValueType = typename InputViewType::value_type; + ValueType _offset; OutputViewType _output; - inputViewType _input; + InputViewType _input; KOKKOS_INLINE_FUNCTION - F_checkPointInclusion(const double offset_, + F_checkPointInclusion(const ValueType offset_, OutputViewType output_, - inputViewType input_) + InputViewType input_) : _offset(offset_), _output(output_), _input(input_) {} @@ -92,9 +89,8 @@ namespace Intrepid2 { void operator()(const ordinal_type i) const { const auto in = Kokkos::subview(_input,i,Kokkos::ALL()); - for (int k=0;k<3;++k) in(k)+=_offset; - const auto check = cellTopologyTagType::checkPointInclusion(in, 0.0); - + for (int k=0;k<3;++k) in(k)+=_offset; + const auto check = PointInclusion::check(in, threshold()); _output(i) = check; } }; @@ -102,67 +98,67 @@ namespace Intrepid2 { // We provide maps belonging to the functional spaces generated by the finite element basis functions associated to the each topology - template + template struct MapPoints; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int /*comp*/) {const auto& x = coords[0]; return 2.0+3*x;} + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int /*comp*/) {const auto& x = coords[0]; return 2.0+3*x;} }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int /*comp*/) {const auto& x = coords[0]; return 2.0+3*x-0.1*x*x;} + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int /*comp*/) {const auto& x = coords[0]; return 2.0+3*x-0.1*x*x;} }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { const auto& x = coords[0]; const auto& y = coords[1]; return (comp == 0) ? 2.0+3*x+2*y : -2.0+2*x+5*y; } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { - MapPoints::key> linearMap; + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { + MapPoints::key, ValueType> linearMap; const auto& x = coords[0]; const auto& y = coords[1]; return linearMap(coords, comp) -0.1*x*x+0.05*x*y+0.2*y*y; } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { const auto& x = coords[0]; const auto& y = coords[1]; return (comp == 0) ? 2.0+3*x+2*y+0.1*x*y : -2.0+2*x+5*y-0.1*x*y; } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { - MapPoints::key> bilinearMap; + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { + MapPoints::key, ValueType> bilinearMap; const auto& x = coords[0]; const auto& y = coords[1]; return bilinearMap(coords, comp)-0.1*x*x -0.2*y*y +0.05*x*y*(x-y); } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { - MapPoints::key> bilinearMap; + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { + MapPoints::key, ValueType> bilinearMap; const auto& x = coords[0]; const auto& y = coords[1]; return bilinearMap(coords, comp)-0.1*x*x -0.2*y*y +0.05*x*y*(x-y+x*y); } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { const auto& x = coords[0]; const auto& y = coords[1]; const auto& z = coords[2]; return (comp == 0) ? 2.0+3*x+2*y+4*z : (comp == 1) ? -2.0+2*x+5*y+4*z : @@ -170,39 +166,39 @@ namespace Intrepid2 { } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { - MapPoints::key> linearMap; + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { + MapPoints::key, ValueType> linearMap; const auto& x = coords[0]; const auto& y = coords[1]; const auto& z = coords[2]; return linearMap(coords, comp)-0.1*x*x -0.2*y*y +0.3*z*z +0.05*x*y +0.07*x*z -0.06*y*z; //for simplicity we choose the same higher-order terms for all components } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { - MapPoints::key> linearMap; - const double eps = epsilon(); + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { + MapPoints::key, ValueType> linearMap; + const ValueType eps = epsilon(); const auto& x = coords[0]; const auto& y = coords[1]; const auto z = (1.0 - coords[2] < eps) ? 1.0-eps : coords[2]; return linearMap(coords, comp)-0.1*x*y/(1.0-z); //for simplicity we choose the same higher-order terms for all components } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { - MapPoints::key> quadraticMap; - const double eps = epsilon(); + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { + MapPoints::key, ValueType> quadraticMap; + const ValueType eps = epsilon(); const auto& x = coords[0]; const auto& y = coords[1]; const auto z = (1.0 - coords[2] < eps) ? 1.0-eps : coords[2]; return quadraticMap(coords, comp)-0.1*x*y/(1.0-z)*(1.0-x+y); //for simplicity we choose the higher-order terms for all components } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { const auto& x = coords[0]; const auto& y = coords[1]; const auto& z = coords[2]; return (comp == 0) ? 2.0+3*x+2*y+4*z +0.07*x*z -0.06*y*z : (comp == 1) ? -2.0+2*x+5*y+4*z +0.07*x*z -0.06*y*z : @@ -210,27 +206,27 @@ namespace Intrepid2 { } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { - MapPoints::key> wedge6Map; + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { + MapPoints::key, ValueType> wedge6Map; const auto& x = coords[0]; const auto& y = coords[1]; const auto& z = coords[2]; return wedge6Map(coords, comp)-0.05*(x*x - y*y + z*z + x*y + x*y*z + y*z*z + x*x*z + y*y*z + x*z*z); //for simplicity we choose the same higher-order terms for all components } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { - MapPoints::key> wedge15Map; + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { + MapPoints::key, ValueType> wedge15Map; const auto& x = coords[0]; const auto& y = coords[1]; const auto& z = coords[2]; return wedge15Map(coords, comp)+0.04*(x*x*z*z - x*y*z*z + y*y*z*z); //for simplicity we choose the same higher-order terms for all components } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { const auto& x = coords[0]; const auto& y = coords[1]; const auto& z = coords[2]; return (comp == 0) ? 2.0+3*x+2*y+4*z +0.05*x*y +0.07*x*z -0.06*y*z + 0.05*x*y*z : (comp == 1) ? -2.0+2*x+5*y+4*z +0.05*x*y +0.07*x*z -0.06*y*z + 0.05*x*y*z : @@ -238,55 +234,51 @@ namespace Intrepid2 { } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { - MapPoints::key> trilinearMap; + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { + MapPoints::key, ValueType> trilinearMap; const auto& x = coords[0]; const auto& y = coords[1]; const auto& z = coords[2]; return trilinearMap(coords, comp)-0.05*(x*x - y*y + z*z + x*y*z*(x-y-z)+ x*x*y + y*z*z + x*x*z + y*y*z+ x*y*y + x*z*z); //for simplicity we choose the same higher-order terms for all components } }; - template <> - struct MapPoints::key> { - double operator()(const double* coords, int comp) { - MapPoints::key> hex20Map; + template + struct MapPoints::key, ValueType> { + ValueType operator()(const ValueType* coords, int comp) { + MapPoints::key, ValueType> hex20Map; const auto& x = coords[0]; const auto& y = coords[1]; const auto& z = coords[2]; return hex20Map(coords, comp)+0.07*(x*x*y*y + x*x*z*z + y*y*z*z + x*x*y*y*z*z + x*y*z*(x*y-x*z+y*z)); //for simplicity we choose the same higher-order terms for all components } }; - - - double mapLine2(const double* coords, int /*comp*/) {const auto& x = coords[0]; return 2.0+3*x;} - // This macro computes the reference cell nodes and maps them to the physical space (according to functional space associated to the topology). // It also maps four input points to the same physical sapce. #define INTREPID2_COMPUTE_POINTS_AND_CELL_NODES_IN_PHYS_SPACE(inCell, host_points, physPoints, physNodes, shtopo, cellTools) \ { \ - MapPoints mapPoints; \ - Intrepid2::HostBasisPtr basisPtr; \ + MapPoints mapPoints; \ + Intrepid2::HostBasisPtr basisPtr; \ using HostDeviceType = Kokkos::HostSpace::device_type; \ shards::CellTopology topo( shards::getCellTopologyData()); \ switch (topo.getKey()) { \ - case shards::Line<2>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_LINE_C1_FEM ()); break; \ - case shards::Line<3>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_LINE_C2_FEM ()); break; \ - case shards::Triangle<3>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_TRI_C1_FEM ()); break; \ - case shards::Quadrilateral<4>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_QUAD_C1_FEM ()); break; \ - case shards::Tetrahedron<4>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_TET_C1_FEM ()); break; \ - case shards::Hexahedron<8>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_HEX_C1_FEM ()); break; \ - case shards::Wedge<6>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_WEDGE_C1_FEM ()); break; \ - case shards::Pyramid<5>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_PYR_C1_FEM ()); break; \ - case shards::Triangle<6>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_TRI_C2_FEM ()); break; \ - case shards::Quadrilateral<8>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_QUAD_I2_FEM ()); break; \ - case shards::Quadrilateral<9>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_QUAD_C2_FEM ()); break; \ - case shards::Tetrahedron<10>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_TET_C2_FEM ()); break; \ - case shards::Tetrahedron<11>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_TET_COMP12_FEM()); break; \ - case shards::Hexahedron<20>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_HEX_I2_FEM ()); break; \ - case shards::Hexahedron<27>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_HEX_C2_FEM ()); break; \ - case shards::Wedge<15>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_WEDGE_I2_FEM ()); break; \ - case shards::Wedge<18>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_WEDGE_C2_FEM ()); break; \ - case shards::Pyramid<13>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_PYR_I2_FEM ()); break; \ + case shards::Line<2>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_LINE_C1_FEM ()); break; \ + case shards::Line<3>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_LINE_C2_FEM ()); break; \ + case shards::Triangle<3>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_TRI_C1_FEM ()); break; \ + case shards::Quadrilateral<4>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_QUAD_C1_FEM ()); break; \ + case shards::Tetrahedron<4>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_TET_C1_FEM ()); break; \ + case shards::Hexahedron<8>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_HEX_C1_FEM ()); break; \ + case shards::Wedge<6>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_WEDGE_C1_FEM ()); break; \ + case shards::Pyramid<5>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_PYR_C1_FEM ()); break; \ + case shards::Triangle<6>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_TRI_C2_FEM ()); break; \ + case shards::Quadrilateral<8>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_QUAD_I2_FEM ()); break; \ + case shards::Quadrilateral<9>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_QUAD_C2_FEM ()); break; \ + case shards::Tetrahedron<10>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_TET_C2_FEM ()); break; \ + case shards::Tetrahedron<11>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_TET_COMP12_FEM()); break; \ + case shards::Hexahedron<20>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_HEX_I2_FEM ()); break; \ + case shards::Hexahedron<27>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_HEX_C2_FEM ()); break; \ + case shards::Wedge<15>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_WEDGE_I2_FEM ()); break; \ + case shards::Wedge<18>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_WEDGE_C2_FEM ()); break; \ + case shards::Pyramid<13>::key: basisPtr = Teuchos::rcp(new Basis_HGRAD_PYR_I2_FEM ()); break; \ default: {}; \ } \ \ @@ -294,10 +286,10 @@ namespace Intrepid2 { auto refNodes_h = Kokkos::DynRankView("refNodes",basisPtr->getCardinality(),dim); \ basisPtr->getDofCoords(refNodes_h); \ physNodes = Kokkos::DynRankView("physNodes",1,basisPtr->getCardinality(),dim); \ - physPoints = Kokkos::DynRankView("physPoints", host_points.extent(0),dim); \ + physPoints = Kokkos::DynRankView("physPoints", host_points.extent(0),dim); \ auto physNodes_h = Kokkos::create_mirror_view(physNodes); \ auto physPoints_h = Kokkos::create_mirror_view(physPoints); \ - double coords[3]={}; \ + ValueType coords[3]={}; \ for(size_t i=0; i< refNodes_h.extent(0);++i) { \ for (size_t d=0; d()); \ assert(ct::hasReferenceCell(topo)); \ std::string label = "reference"; \ @@ -379,7 +372,7 @@ namespace Intrepid2 { << "| |\n" << "===============================================================================\n"; - const ValueType tol = tolerence()*100.0; + const ValueType tol = tolerence()*100.0; int errorFlag = 0; @@ -392,30 +385,30 @@ namespace Intrepid2 { << "===============================================================================\n\n"; { - double offset = 0.0; - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Line<>, Impl::Line<2>); + ValueType offset = 0.0; + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Line<2>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Triangle<>, Impl::Triangle<3>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Quadrilateral<>, Impl::Quadrilateral<4>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Triangle<3>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Quadrilateral<4>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Tetrahedron<>, Impl::Tetrahedron<4>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Hexahedron<>, Impl::Hexahedron<8>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Tetrahedron<4>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Hexahedron<8>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Pyramid<>, Impl::Pyramid<5>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Wedge<>, Impl::Wedge<6>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Pyramid<5>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, true, shards::Wedge<6>); } { - double offset = 3.0; - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Line<>, Impl::Line<2>); + ValueType offset = 3.0; + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Line<2>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Triangle<>, Impl::Triangle<3>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Quadrilateral<>, Impl::Quadrilateral<4>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Triangle<3>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Quadrilateral<4>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Tetrahedron<>, Impl::Tetrahedron<4>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Hexahedron<>, Impl::Hexahedron<8>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Tetrahedron<4>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Hexahedron<8>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Pyramid<>, Impl::Pyramid<5>); - INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Wedge<>, Impl::Wedge<6>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Pyramid<5>); + INTREPID2_TEST_CHECK_POINT_INCLUSION(offset, false, shards::Wedge<6>); } @@ -433,7 +426,7 @@ namespace Intrepid2 { *outStream << "\n" << "===============================================================================\n" - << "| Test 2: test Point wise inclusion points\n" + << "| Test 2: test Point wise inclusion \n" << "===============================================================================\n\n"; @@ -453,7 +446,7 @@ namespace Intrepid2 { auto pts1d_h = Kokkos::create_mirror_view(pts1d); auto pts2d_h = Kokkos::create_mirror_view(pts2d); auto pts3d_h = Kokkos::create_mirror_view(pts3d); - double eps = 1e-4; + ValueType eps = 1e-4; using ct = Intrepid2::CellTools; // line topologies