49#ifndef Intrepid2_TestUtils_h
50#define Intrepid2_TestUtils_h
52#include "Kokkos_Core.hpp"
53#include "Kokkos_DynRankView.hpp"
63#include "Teuchos_UnitTestHarness.hpp"
71#if defined(KOKKOS_ENABLE_CUDA)
73#elif defined(KOKKOS_ENABLE_HIP)
74 using DefaultTestDeviceType = Kokkos::Device<Kokkos::Experimental::HIP,Kokkos::Experimental::HIPSpace>;
80 template <
class Scalar>
81 const typename Teuchos::ScalarTraits<Scalar>::magnitudeType
84 using ST = Teuchos::ScalarTraits<Scalar>;
85 return ST::magnitude(Teuchos::RelErrSmallNumber<ST::hasMachineParameters,Scalar>::smallNumber());
89 template <
class Scalar1,
class Scalar2>
90 KOKKOS_INLINE_FUNCTION
92 relErrMeetsTol(
const Scalar1 &s1,
const Scalar2 &s2,
const typename Teuchos::ScalarTraits<
typename std::common_type<Scalar1,Scalar2>::type >::magnitudeType &smallNumber,
const double &tol )
94 using Scalar =
typename std::common_type<Scalar1,Scalar2>::type;
95 const Scalar s1Abs = fabs(s1);
96 const Scalar s2Abs = fabs(s2);
97 const Scalar maxAbs = (s1Abs > s2Abs) ? s1Abs : s2Abs;
98 const Scalar relErr = fabs( s1 - s2 ) / (
smallNumber + maxAbs );
102 template <
class Scalar1,
class Scalar2>
103 KOKKOS_INLINE_FUNCTION
105 errMeetsAbsAndRelTol(
const Scalar1 &s1,
const Scalar2 &s2,
const double &relTol,
const double &absTol )
107 return fabs( s1 - s2 ) <= absTol + fabs(s1) * relTol;
110 static const double TEST_TOLERANCE_TIGHT = 1.e2 * std::numeric_limits<double>::epsilon();
113 template<
typename ScalarType,
typename DeviceType>
114 using ViewType = Kokkos::DynRankView<ScalarType,DeviceType>;
116 template<
typename ScalarType,
typename DeviceType>
117 using FixedRankViewType = Kokkos::View<ScalarType,DeviceType>;
119 template<
typename ScalarType>
120 KOKKOS_INLINE_FUNCTION
bool valuesAreSmall(
const ScalarType &a,
const ScalarType &b,
const double &epsilon)
123 return (abs(a) < epsilon) && (abs(b) < epsilon);
126 inline bool approximatelyEqual(
double a,
double b,
double epsilon)
128 const double larger_magnitude = (std::abs(a) < std::abs(b) ? std::abs(b) : std::abs(a));
129 return std::abs(a - b) <= larger_magnitude * epsilon;
132 inline bool essentiallyEqual(
double a,
double b,
double epsilon)
134 const double smaller_magnitude = (std::abs(a) > std::abs(b) ? std::abs(b) : std::abs(a));
135 return std::abs(a - b) <= smaller_magnitude * epsilon;
139 KOKKOS_INLINE_FUNCTION
double fromZeroOne(
double x_zero_one)
141 return x_zero_one * 2.0 - 1.0;
145 KOKKOS_INLINE_FUNCTION
double toZeroOne(
double x_minus_one_one)
147 return (x_minus_one_one + 1.0) / 2.0;
151 KOKKOS_INLINE_FUNCTION
double fromZeroOne_dx(
double dx_zero_one)
153 return dx_zero_one / 2.0;
157 KOKKOS_INLINE_FUNCTION
double toZeroOne_dx(
double dx_minus_one_one)
159 return dx_minus_one_one * 2.0;
162 template<
class DeviceViewType>
163 typename DeviceViewType::HostMirror getHostCopy(
const DeviceViewType &deviceView)
165 typename DeviceViewType::HostMirror hostView = Kokkos::create_mirror(deviceView);
166 Kokkos::deep_copy(hostView, deviceView);
170 template<
class BasisFamily>
171 inline Teuchos::RCP< Intrepid2::Basis<DefaultTestDeviceType,double,double> > getBasisUsingFamily(shards::CellTopology cellTopo, Intrepid2::EFunctionSpace fs,
172 int polyOrder_x,
int polyOrder_y=-1,
int polyOrder_z = -1)
174 using BasisPtr =
typename BasisFamily::BasisPtr;
177 using namespace Intrepid2;
179 if (cellTopo.getBaseKey() == shards::Line<>::key)
181 basis = getLineBasis<BasisFamily>(fs,polyOrder_x);
183 else if (cellTopo.getBaseKey() == shards::Quadrilateral<>::key)
185 INTREPID2_TEST_FOR_EXCEPTION(polyOrder_y < 0, std::invalid_argument,
"polyOrder_y must be specified");
186 basis = getQuadrilateralBasis<BasisFamily>(fs,polyOrder_x,polyOrder_y);
188 else if (cellTopo.getBaseKey() == shards::Triangle<>::key)
190 basis = getTriangleBasis<BasisFamily>(fs,polyOrder_x);
192 else if (cellTopo.getBaseKey() == shards::Hexahedron<>::key)
194 INTREPID2_TEST_FOR_EXCEPTION(polyOrder_y < 0, std::invalid_argument,
"polyOrder_y must be specified");
195 INTREPID2_TEST_FOR_EXCEPTION(polyOrder_z < 0, std::invalid_argument,
"polyOrder_z must be specified");
196 basis = getHexahedronBasis<BasisFamily>(fs,polyOrder_x,polyOrder_y,polyOrder_z);
198 else if (cellTopo.getBaseKey() == shards::Tetrahedron<>::key)
200 basis = getTetrahedronBasis<BasisFamily>(fs, polyOrder_x);
202 else if (cellTopo.getBaseKey() == shards::Wedge<>::key)
204 INTREPID2_TEST_FOR_EXCEPTION(polyOrder_y < 0, std::invalid_argument,
"polyOrder_y must be specified");
205 basis = getWedgeBasis<BasisFamily>(fs,polyOrder_x,polyOrder_y);
209 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported cell topology");
214 template<
bool defineVertexFunctions>
215 inline Teuchos::RCP< Intrepid2::Basis<DefaultTestDeviceType,double,double> > getHierarchicalBasis(shards::CellTopology cellTopo, Intrepid2::EFunctionSpace fs,
216 int polyOrder_x,
int polyOrder_y=-1,
int polyOrder_z = -1)
219 using Scalar = double;
220 using namespace Intrepid2;
229 return getBasisUsingFamily<BasisFamily>(cellTopo, fs, polyOrder_x, polyOrder_y, polyOrder_z);
232 template<
typename ValueType,
typename DeviceType,
class ... DimArgs>
233 inline ViewType<ValueType,DeviceType> getView(
const std::string &label, DimArgs... dims)
235 const bool allocateFadStorage = !std::is_pod<ValueType>::value;
236 if (!allocateFadStorage)
238 return ViewType<ValueType,DeviceType>(label,dims...);
242 return ViewType<ValueType,DeviceType>(label,dims...,MAX_FAD_DERIVATIVES_FOR_TESTS+1);
247 template<
typename ValueType,
class ... DimArgs>
248 inline FixedRankViewType<
typename RankExpander<ValueType,
sizeof...(DimArgs) >::value_type,
DefaultTestDeviceType > getFixedRankView(
const std::string &label, DimArgs... dims)
250 const bool allocateFadStorage = !std::is_pod<ValueType>::value;
251 using value_type =
typename RankExpander<ValueType,
sizeof...(dims) >::value_type;
252 if (!allocateFadStorage)
254 return FixedRankViewType<value_type,DefaultTestDeviceType>(label,dims...);
258 return FixedRankViewType<value_type,DefaultTestDeviceType>(label,dims...,MAX_FAD_DERIVATIVES_FOR_TESTS+1);
268 template <
typename Po
intValueType,
typename DeviceType>
269 inline ViewType<PointValueType,DeviceType>
getInputPointsView(shards::CellTopology &cellTopo,
int numPoints_1D)
271 if (cellTopo.getBaseKey() == shards::Wedge<>::key)
273 shards::CellTopology lineTopo = shards::CellTopology(shards::getCellTopologyData<shards::Line<> >() );
274 shards::CellTopology triTopo = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >() );
276 const ordinal_type order = numPoints_1D - 1;
279 ordinal_type numPoints = numPoints_tri * numPoints_line;
280 ordinal_type spaceDim = cellTopo.getDimension();
282 ViewType<PointValueType,DeviceType> inputPointsTri = getView<PointValueType,DeviceType>(
"input points",numPoints_tri, 2);
283 ViewType<PointValueType,DeviceType> inputPointsLine = getView<PointValueType,DeviceType>(
"input points",numPoints_line,1);
287 ViewType<PointValueType,DeviceType> inputPoints = getView<PointValueType,DeviceType>(
"input points",numPoints,spaceDim);
289 using ExecutionSpace =
typename ViewType<PointValueType,DeviceType>::execution_space;
291 Kokkos::RangePolicy < ExecutionSpace > policy(0,numPoints_tri);
292 Kokkos::parallel_for( policy,
293 KOKKOS_LAMBDA (
const ordinal_type &triPointOrdinal )
295 ordinal_type pointOrdinal = triPointOrdinal * numPoints_line;
296 for (ordinal_type linePointOrdinal=0; linePointOrdinal<numPoints_line; linePointOrdinal++)
298 inputPoints(pointOrdinal,0) = inputPointsTri( triPointOrdinal,0);
299 inputPoints(pointOrdinal,1) = inputPointsTri( triPointOrdinal,1);
300 inputPoints(pointOrdinal,2) = inputPointsLine(linePointOrdinal,0);
310 const ordinal_type order = numPoints_1D - 1;
312 ordinal_type spaceDim = cellTopo.getDimension();
314 ViewType<PointValueType,DeviceType> inputPoints = getView<PointValueType,DeviceType>(
"input points",numPoints,spaceDim);
321 template<
typename OutputValueType,
typename DeviceType>
322 inline ViewType<OutputValueType,DeviceType> getOutputView(Intrepid2::EFunctionSpace fs, Intrepid2::EOperator op,
int basisCardinality,
int numPoints,
int spaceDim)
325 case Intrepid2::FUNCTION_SPACE_HGRAD:
327 case Intrepid2::OPERATOR_VALUE:
328 return getView<OutputValueType,DeviceType>(
"H^1 value output",basisCardinality,numPoints);
329 case Intrepid2::OPERATOR_GRAD:
330 return getView<OutputValueType,DeviceType>(
"H^1 derivative output",basisCardinality,numPoints,spaceDim);
331 case Intrepid2::OPERATOR_D1:
332 case Intrepid2::OPERATOR_D2:
333 case Intrepid2::OPERATOR_D3:
334 case Intrepid2::OPERATOR_D4:
335 case Intrepid2::OPERATOR_D5:
336 case Intrepid2::OPERATOR_D6:
337 case Intrepid2::OPERATOR_D7:
338 case Intrepid2::OPERATOR_D8:
339 case Intrepid2::OPERATOR_D9:
340 case Intrepid2::OPERATOR_D10:
342 const auto dkcard = getDkCardinality(op, spaceDim);
343 return getView<OutputValueType,DeviceType>(
"H^1 derivative output",basisCardinality,numPoints,dkcard);
346 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported op/fs combination");
348 case Intrepid2::FUNCTION_SPACE_HCURL:
350 case Intrepid2::OPERATOR_VALUE:
351 return getView<OutputValueType,DeviceType>(
"H(curl) value output",basisCardinality,numPoints,spaceDim);
352 case Intrepid2::OPERATOR_CURL:
354 return getView<OutputValueType,DeviceType>(
"H(curl) derivative output",basisCardinality,numPoints);
355 else if (spaceDim == 3)
356 return getView<OutputValueType,DeviceType>(
"H(curl) derivative output",basisCardinality,numPoints,spaceDim);
358 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported op/fs combination");
360 case Intrepid2::FUNCTION_SPACE_HDIV:
362 case Intrepid2::OPERATOR_VALUE:
363 return getView<OutputValueType,DeviceType>(
"H(div) value output",basisCardinality,numPoints,spaceDim);
364 case Intrepid2::OPERATOR_DIV:
365 return getView<OutputValueType,DeviceType>(
"H(div) derivative output",basisCardinality,numPoints);
367 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported op/fs combination");
370 case Intrepid2::FUNCTION_SPACE_HVOL:
372 case Intrepid2::OPERATOR_VALUE:
373 return getView<OutputValueType,DeviceType>(
"H(vol) value output",basisCardinality,numPoints);
374 case Intrepid2::OPERATOR_D1:
375 case Intrepid2::OPERATOR_D2:
376 case Intrepid2::OPERATOR_D3:
377 case Intrepid2::OPERATOR_D4:
378 case Intrepid2::OPERATOR_D5:
379 case Intrepid2::OPERATOR_D6:
380 case Intrepid2::OPERATOR_D7:
381 case Intrepid2::OPERATOR_D8:
382 case Intrepid2::OPERATOR_D9:
383 case Intrepid2::OPERATOR_D10:
386 return getView<OutputValueType,DeviceType>(
"H(vol) derivative output",basisCardinality,numPoints,dkcard);
389 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported op/fs combination");
392 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported op/fs combination");
398 inline std::vector< std::vector<int> > getBasisTestCasesUpToDegree(
int spaceDim,
int minDegree,
int polyOrder_x,
int polyOrder_y=-1,
int polyOrder_z=-1)
400 std::vector<int> degrees(spaceDim);
401 degrees[0] = polyOrder_x;
402 if (spaceDim > 1) degrees[1] = polyOrder_y;
403 if (spaceDim > 2) degrees[2] = polyOrder_z;
405 int numCases = degrees[0];
406 for (
unsigned d=1; d<degrees.size(); d++)
408 INTREPID2_TEST_FOR_EXCEPTION(degrees[d] < minDegree, std::invalid_argument,
"Unsupported degree/minDegree combination");
409 numCases = numCases * (degrees[d] + 1 - minDegree);
411 std::vector< std::vector<int> > subBasisDegreeTestCases(numCases);
412 for (
int caseOrdinal=0; caseOrdinal<numCases; caseOrdinal++)
414 std::vector<int> subBasisDegrees(degrees.size());
415 int caseRemainder = caseOrdinal;
416 for (
int d=degrees.size()-1; d>=0; d--)
418 int subBasisDegree = caseRemainder % (degrees[d] + 1 - minDegree);
419 caseRemainder = caseRemainder / (degrees[d] + 1 - minDegree);
420 subBasisDegrees[d] = subBasisDegree + minDegree;
422 subBasisDegreeTestCases[caseOrdinal] = subBasisDegrees;
424 return subBasisDegreeTestCases;
428 template<
class Functor,
class Scalar,
int rank>
431 INTREPID2_TEST_FOR_EXCEPTION(rank !=
getFunctorRank(deviceFunctor), std::invalid_argument,
"functor rank must match the template argument");
434 ViewType<Scalar,DeviceType> view;
435 const std::string label =
"functor copy";
436 const auto &f = deviceFunctor;
440 view = getView<Scalar,DeviceType>(label);
443 view = getView<Scalar,DeviceType>(label, f.extent_int(0));
446 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1));
449 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2));
452 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2), f.extent_int(3));
455 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2), f.extent_int(3), f.extent_int(4));
458 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2), f.extent_int(3), f.extent_int(4), f.extent_int(5));
461 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2), f.extent_int(3), f.extent_int(4), f.extent_int(5), f.extent_int(6));
464 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported functor rank");
467 int entryCount = view.size();
469 using ExecutionSpace =
typename ViewType<Scalar,DeviceType>::execution_space;
474 Kokkos::RangePolicy < ExecutionSpace > policy(0,entryCount);
475 Kokkos::parallel_for( policy,
476 KOKKOS_LAMBDA (
const int &enumerationIndex )
478 ViewIteratorScalar vi(view);
479 FunctorIteratorScalar fi(f);
481 vi.setEnumerationIndex(enumerationIndex);
482 fi.setEnumerationIndex(enumerationIndex);
488 auto hostView = Kokkos::create_mirror_view(view);
489 Kokkos::deep_copy(hostView, view);
493 template<
class FunctorType,
typename Scalar,
int rank>
494 void printFunctor(
const FunctorType &functor, std::ostream &out,
const std::string &functorName =
"")
496 auto functorHostCopy = copyFunctorToHostView<FunctorType, Scalar, rank>(functor);
498 std::string name = (functorName ==
"") ?
"Functor" : functorName;
500 out <<
"\n******** " << name <<
" (rank " << rank <<
") ********\n";
501 out <<
"dimensions: (";
502 for (
int r=0; r<rank; r++)
504 out << functor.extent_int(r);
505 if (r<rank-1) out <<
",";
509 ViewIterator<
decltype(functorHostCopy),Scalar> vi(functorHostCopy);
511 bool moreEntries =
true;
514 Scalar value = vi.get();
516 auto location = vi.getLocation();
517 out << functorName <<
"(";
518 for (ordinal_type i=0; i<rank; i++)
526 out <<
") " << value << std::endl;
528 moreEntries = (vi.increment() != -1);
533 template<
class FunctorType>
534 void printFunctor1(
const FunctorType &functor, std::ostream &out,
const std::string &functorName =
"")
536 using Scalar =
typename std::remove_reference<
decltype(functor(0))>::type;
537 printFunctor<FunctorType, Scalar, 1>(functor, out, functorName);
540 template<
class FunctorType>
541 void printFunctor2(
const FunctorType &functor, std::ostream &out,
const std::string &functorName =
"")
543 using Scalar =
typename std::remove_const<
typename std::remove_reference<
decltype(functor(0,0))>::type>::type;
544 printFunctor<FunctorType, Scalar, 2>(functor, out, functorName);
547 template<
class FunctorType>
548 void printFunctor3(
const FunctorType &functor, std::ostream &out,
const std::string &functorName =
"")
550 using Scalar =
typename std::remove_const<
typename std::remove_reference<
decltype(functor(0,0,0))>::type>::type;
551 printFunctor<FunctorType, Scalar, 3>(functor, out, functorName);
554 template<
class FunctorType>
555 void printFunctor4(
const FunctorType &functor, std::ostream &out,
const std::string &functorName =
"")
557 using Scalar =
typename std::remove_const<
typename std::remove_reference<
decltype(functor(0,0,0,0))>::type>::type;
558 printFunctor<FunctorType, Scalar, 4>(functor, out, functorName);
561 template<
class FunctorType>
562 void printFunctor5(
const FunctorType &functor, std::ostream &out,
const std::string &functorName =
"")
564 using Scalar =
typename std::remove_const<
typename std::remove_reference<
decltype(functor(0,0,0,0,0))>::type>::type;
565 printFunctor<FunctorType, Scalar, 5>(functor, out, functorName);
568 template<
class FunctorType>
569 void printFunctor6(
const FunctorType &functor, std::ostream &out,
const std::string &functorName =
"")
571 using Scalar =
typename std::remove_const<
typename std::remove_reference<
decltype(functor(0,0,0,0,0,0))>::type>::type;
572 printFunctor<FunctorType, Scalar, 6>(functor, out, functorName);
575 template<
class FunctorType>
576 void printFunctor7(
const FunctorType &functor, std::ostream &out,
const std::string &functorName =
"")
578 using Scalar =
typename std::remove_const<
typename std::remove_reference<
decltype(functor(0,0,0,0,0,0,0))>::type>::type;
579 printFunctor<FunctorType, Scalar, 7>(functor, out, functorName);
583 void printView(
const View &view, std::ostream &out,
const std::string &viewName =
"")
585 using Scalar =
typename View::value_type;
586 using HostView =
typename View::HostMirror;
589 auto hostView = getHostCopy(view);
591 HostViewIteratorScalar vi(hostView);
593 bool moreEntries = (vi.nextIncrementRank() != -1);
596 Scalar value = vi.get();
598 auto location = vi.getLocation();
599 out << viewName <<
"(";
608 out <<
") " << value << std::endl;
610 moreEntries = (vi.increment() != -1);
614 template <
class FunctorType1,
class FunctorType2,
int rank,
typename Scalar=
typename FunctorType1::value_type,
class ExecutionSpace =
typename FunctorType1::execution_space>
616 testFloatingEquality(
const FunctorType1 &functor1,
const FunctorType2 &functor2,
double relTol,
double absTol, Teuchos::FancyOStream &out,
bool &success,
617 std::string functor1Name =
"Functor 1", std::string functor2Name =
"Functor 2")
619 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"testFloatingEquality() called on FunctorType1 or FunctorType2 that does not support the specified rank.\n");
623 template <
class FunctorType1,
class FunctorType2,
int rank,
typename Scalar=
typename FunctorType1::value_type,
class ExecutionSpace =
typename FunctorType1::execution_space>
625 testFloatingEquality(
const FunctorType1 &functor1,
const FunctorType2 &functor2,
double relTol,
double absTol, Teuchos::FancyOStream &out,
bool &success,
626 std::string functor1Name =
"Functor 1", std::string functor2Name =
"Functor 2")
635 TEUCHOS_TEST_FOR_EXCEPTION(
getFunctorRank(functor1) != rank, std::invalid_argument,
"functor1 and functor2 must agree in rank");
636 TEUCHOS_TEST_FOR_EXCEPTION(
getFunctorRank(functor2) != rank, std::invalid_argument,
"functor1 and functor2 must agree in rank");
641 TEUCHOS_TEST_FOR_EXCEPTION(functor1.extent_int(i) != functor2.extent_int(i), std::invalid_argument,
642 "functor1 and functor2 must agree in size in each dimension; functor1 has extent "
643 + std::to_string(functor1.extent_int(i)) +
" in dimension " + std::to_string(i)
644 +
"; functor2 has extent " + std::to_string(functor2.extent_int(i)) );
645 entryCount *= functor1.extent_int(i);
647 if (entryCount == 0)
return;
649 ViewType<bool,ExecutionSpace> valuesMatch = getView<bool,ExecutionSpace>(
"valuesMatch", entryCount);
651 Kokkos::RangePolicy < ExecutionSpace > policy(0,entryCount);
652 Kokkos::parallel_for( policy,
653 KOKKOS_LAMBDA (
const int &enumerationIndex )
655 Functor1IteratorScalar vi1(functor1);
656 Functor2IteratorScalar vi2(functor2);
658 vi1.setEnumerationIndex(enumerationIndex);
659 vi2.setEnumerationIndex(enumerationIndex);
661 const Scalar & value1 = vi1.get();
662 const Scalar & value2 = vi2.get();
664 bool errMeetsTol = errMeetsAbsAndRelTol(value1, value2, relTol, absTol);
665 valuesMatch(enumerationIndex) = errMeetsTol;
669 int failureCount = 0;
670 Kokkos::RangePolicy<ExecutionSpace > reducePolicy(0, entryCount);
671 Kokkos::parallel_reduce( reducePolicy,
672 KOKKOS_LAMBDA(
const int &enumerationIndex,
int &reducedValue )
674 reducedValue += valuesMatch(enumerationIndex) ? 0 : 1;
677 const bool allValuesMatch = (failureCount == 0);
678 success = success && allValuesMatch;
683 auto functor1CopyHost = copyFunctorToHostView<FunctorType1,Scalar,rank>(functor1);
684 auto functor2CopyHost = copyFunctorToHostView<FunctorType2,Scalar,rank>(functor2);
686 auto valuesMatchHost = getHostCopy(valuesMatch);
688 FunctorIterator<
decltype(functor1CopyHost),Scalar,rank> vi1(functor1CopyHost);
689 FunctorIterator<
decltype(functor2CopyHost),Scalar,rank> vi2(functor2CopyHost);
692 bool moreEntries =
true;
695 bool errMeetsTol = viMatch.get();
699 const Scalar value1 = vi1.get();
700 const Scalar value2 = vi2.get();
703 auto location = vi1.getLocation();
704 out <<
"At location (";
713 out <<
") " << functor1Name <<
" value != " << functor2Name <<
" value (";
714 out << value1 <<
" != " << value2 <<
"); diff is " << std::abs(value1-value2) << std::endl;
717 moreEntries = (vi1.increment() != -1);
718 moreEntries = moreEntries && (vi2.increment() != -1);
719 moreEntries = moreEntries && (viMatch.increment() != -1);
724 template <
class ViewType,
class FunctorType>
725 void testFloatingEquality1(
const ViewType &view,
const FunctorType &functor,
double relTol,
double absTol, Teuchos::FancyOStream &out,
bool &success,
726 std::string view1Name =
"View", std::string view2Name =
"Functor")
728 testFloatingEquality<ViewType, FunctorType, 1>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
731 template <
class ViewType,
class FunctorType>
732 void testFloatingEquality2(
const ViewType &view,
const FunctorType &functor,
double relTol,
double absTol, Teuchos::FancyOStream &out,
bool &success,
733 std::string view1Name =
"View", std::string view2Name =
"Functor")
735 testFloatingEquality<ViewType, FunctorType, 2>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
738 template <
class ViewType,
class FunctorType>
739 void testFloatingEquality3(
const ViewType &view,
const FunctorType &functor,
double relTol,
double absTol, Teuchos::FancyOStream &out,
bool &success,
740 std::string view1Name =
"View", std::string view2Name =
"Functor")
742 testFloatingEquality<ViewType, FunctorType, 3>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
745 template <
class ViewType,
class FunctorType>
746 void testFloatingEquality4(
const ViewType &view,
const FunctorType &functor,
double relTol,
double absTol, Teuchos::FancyOStream &out,
bool &success,
747 std::string view1Name =
"View", std::string view2Name =
"Functor")
749 testFloatingEquality<ViewType, FunctorType, 4>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
752 template <
class ViewType,
class FunctorType>
753 void testFloatingEquality5(
const ViewType &view,
const FunctorType &functor,
double relTol,
double absTol, Teuchos::FancyOStream &out,
bool &success,
754 std::string view1Name =
"View", std::string view2Name =
"Functor")
756 testFloatingEquality<ViewType, FunctorType, 5>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
759 template <
class ViewType,
class FunctorType>
760 void testFloatingEquality6(
const ViewType &view,
const FunctorType &functor,
double relTol,
double absTol, Teuchos::FancyOStream &out,
bool &success,
761 std::string view1Name =
"View", std::string view2Name =
"Functor")
763 testFloatingEquality<ViewType, FunctorType, 6>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
766 template <
class ViewType,
class FunctorType>
767 void testFloatingEquality7(
const ViewType &view,
const FunctorType &functor,
double relTol,
double absTol, Teuchos::FancyOStream &out,
bool &success,
768 std::string view1Name =
"View", std::string view2Name =
"Functor")
770 testFloatingEquality<ViewType, FunctorType, 7>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
773 template <
class ViewType1,
class ViewType2>
774 void testViewFloatingEquality(
const ViewType1 &view1,
const ViewType2 &view2,
double relTol,
double absTol, Teuchos::FancyOStream &out,
bool &success,
775 std::string view1Name =
"View 1", std::string view2Name =
"View 2")
778 TEUCHOS_TEST_FOR_EXCEPTION(view1.rank() != view2.rank(), std::invalid_argument,
"views must agree in rank");
779 for (
unsigned i=0; i<view1.rank(); i++)
781 TEUCHOS_TEST_FOR_EXCEPTION(view1.extent_int(i) != view2.extent_int(i), std::invalid_argument,
"views must agree in size in each dimension");
784 if (view1.size() == 0)
return;
786 const int rank = view1.rank();
789 case 0: testFloatingEquality<ViewType1, ViewType2, 0>(view1, view2, relTol, absTol, out, success, view1Name, view2Name);
break;
790 case 1: testFloatingEquality<ViewType1, ViewType2, 1>(view1, view2, relTol, absTol, out, success, view1Name, view2Name);
break;
791 case 2: testFloatingEquality<ViewType1, ViewType2, 2>(view1, view2, relTol, absTol, out, success, view1Name, view2Name);
break;
792 case 3: testFloatingEquality<ViewType1, ViewType2, 3>(view1, view2, relTol, absTol, out, success, view1Name, view2Name);
break;
793 case 4: testFloatingEquality<ViewType1, ViewType2, 4>(view1, view2, relTol, absTol, out, success, view1Name, view2Name);
break;
794 case 5: testFloatingEquality<ViewType1, ViewType2, 5>(view1, view2, relTol, absTol, out, success, view1Name, view2Name);
break;
795 case 6: testFloatingEquality<ViewType1, ViewType2, 6>(view1, view2, relTol, absTol, out, success, view1Name, view2Name);
break;
796 case 7: testFloatingEquality<ViewType1, ViewType2, 7>(view1, view2, relTol, absTol, out, success, view1Name, view2Name);
break;
797 default: INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported rank");
803#ifdef HAVE_INTREPID2_SACADO
804using Sacado_Fad_DFadType = Sacado::Fad::DFad<double>;
805#define INTREPID2_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
807TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( GROUP_NAME, TEST_NAME, double ) \
809TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( GROUP_NAME, TEST_NAME, Sacado_Fad_DFadType ) \
811#define INTREPID2_OUTPUTSCALAR_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
813TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( GROUP_NAME, TEST_NAME, double, double ) \
815TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( GROUP_NAME, TEST_NAME, Sacado_Fad_DFadType, Sacado_Fad_DFadType ) \
818#define INTREPID2_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
820TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( GROUP_NAME, TEST_NAME, double ) \
822#define INTREPID2_OUTPUTSCALAR_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
824TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( GROUP_NAME, TEST_NAME, double, double ) \
Header file for the abstract base class Intrepid2::Basis.
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
KOKKOS_INLINE_FUNCTION ordinal_type getDkCardinality(const EOperator operatorType, const ordinal_type spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
Stateless class representing a family of basis functions, templated on H(vol) and H(grad) on the line...
Defines the Intrepid2::FunctorIterator class, as well as the Intrepid2::functor_returns_ref SFINAE he...
Stateless classes that act as factories for two families of hierarchical bases. HierarchicalBasisFami...
Header file to include all Sacado headers that are required if using Intrepid2 with Sacado types.
ViewType< Scalar, DefaultTestDeviceType >::HostMirror copyFunctorToHostView(const Functor &deviceFunctor)
Copy the values for the specified functor.
ViewType< PointValueType, DeviceType > getInputPointsView(shards::CellTopology &cellTopo, int numPoints_1D)
Returns a DynRankView containing regularly-spaced points on the specified cell topology.
typename Kokkos::DefaultExecutionSpace::device_type DefaultTestDeviceType
Default Kokkos::Device to use for tests; depends on platform.
constexpr int MAX_FAD_DERIVATIVES_FOR_TESTS
Maximum number of derivatives to track for Fad types in tests.
const Teuchos::ScalarTraits< Scalar >::magnitudeType smallNumber()
Use Teuchos small number determination on host; pass this to Intrepid2::relErr() on device.
KOKKOS_INLINE_FUNCTION bool relErrMeetsTol(const Scalar1 &s1, const Scalar2 &s2, const typename Teuchos::ScalarTraits< typename std::common_type< Scalar1, Scalar2 >::type >::magnitudeType &smallNumber, const double &tol)
Adapted from Teuchos::relErr(); for use in code that may be executed on device.
Header function for Intrepid2::Util class and other utility functions.
enable_if_t< has_rank_method< Functor >::value, unsigned > KOKKOS_INLINE_FUNCTION getFunctorRank(const Functor &functor)
A family of basis functions, constructed from H(vol) and H(grad) bases on the line.
essentially, a read-only variant of ViewIterator, for a general functor (extent_int() and rank() supp...
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line.
Basis defining Legendre basis on the line, a polynomial subspace of L^2 (a.k.a. H(vol)) on the line.
A helper class that allows iteration over some part of a Kokkos View, while allowing the calling code...
SFINAE helper to detect whether a type supports a rank-integral-argument operator().
Helper to get Scalar[*+] where the number of *'s matches the given rank.