Intrepid2
Intrepid2_TestUtils.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid2 Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Kyungjoo Kim (kyukim@sandia.gov),
38 // Mauro Perego (mperego@sandia.gov), or
39 // Nate Roberts (nvrober@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
49 #ifndef Intrepid2_TestUtils_h
50 #define Intrepid2_TestUtils_h
51 
52 #include "Kokkos_Core.hpp"
53 #include "Kokkos_DynRankView.hpp"
54 
55 #include "Intrepid2_Basis.hpp"
59 #include "Intrepid2_PointTools.hpp"
60 #include "Intrepid2_Sacado.hpp" // Sacado includes, guarded by the appropriate preprocessor variable
61 #include "Intrepid2_Utils.hpp"
62 
63 #include "Teuchos_UnitTestHarness.hpp"
64 
65 namespace Intrepid2
66 {
68  constexpr int MAX_FAD_DERIVATIVES_FOR_TESTS = 3;
69 
71 #ifdef KOKKOS_ENABLE_CUDA
72  using DefaultTestDeviceType = Kokkos::Device<Kokkos::Cuda,Kokkos::CudaSpace>;
73 #else
74  using DefaultTestDeviceType = typename Kokkos::DefaultExecutionSpace::device_type;
75 #endif
76 
78  template <class Scalar>
79  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType
81  {
82  using ST = Teuchos::ScalarTraits<Scalar>;
83  return ST::magnitude(Teuchos::RelErrSmallNumber<ST::hasMachineParameters,Scalar>::smallNumber());
84  }
85 
87  template <class Scalar1, class Scalar2>
88  KOKKOS_INLINE_FUNCTION
89  bool
90  relErrMeetsTol( const Scalar1 &s1, const Scalar2 &s2, const typename Teuchos::ScalarTraits< typename std::common_type<Scalar1,Scalar2>::type >::magnitudeType &smallNumber, const double &tol )
91  {
92  using Scalar = typename std::common_type<Scalar1,Scalar2>::type;
93  const Scalar s1Abs = fabs(s1);
94  const Scalar s2Abs = fabs(s2);
95  const Scalar maxAbs = (s1Abs > s2Abs) ? s1Abs : s2Abs;
96  const Scalar relErr = fabs( s1 - s2 ) / ( smallNumber + maxAbs );
97  return relErr < tol;
98  }
99 
100  template <class Scalar1, class Scalar2>
101  KOKKOS_INLINE_FUNCTION
102  bool
103  errMeetsAbsAndRelTol( const Scalar1 &s1, const Scalar2 &s2, const double &relTol, const double &absTol )
104  {
105  return fabs( s1 - s2 ) <= absTol + fabs(s1) * relTol;
106  }
107 
108  static const double TEST_TOLERANCE_TIGHT = 1.e2 * std::numeric_limits<double>::epsilon();
109 
110  // we use DynRankView for both input points and values
111  template<typename ScalarType, typename DeviceType>
112  using ViewType = Kokkos::DynRankView<ScalarType,DeviceType>;
113 
114  template<typename ScalarType, typename DeviceType>
115  using FixedRankViewType = Kokkos::View<ScalarType,DeviceType>;
116 
117  template<typename ScalarType>
118  KOKKOS_INLINE_FUNCTION bool valuesAreSmall(const ScalarType &a, const ScalarType &b, const double &epsilon)
119  {
120  using std::abs;
121  return (abs(a) < epsilon) && (abs(b) < epsilon);
122  }
123 
124  inline bool approximatelyEqual(double a, double b, double epsilon)
125  {
126  const double larger_magnitude = (std::abs(a) < std::abs(b) ? std::abs(b) : std::abs(a));
127  return std::abs(a - b) <= larger_magnitude * epsilon;
128  }
129 
130  inline bool essentiallyEqual(double a, double b, double epsilon)
131  {
132  const double smaller_magnitude = (std::abs(a) > std::abs(b) ? std::abs(b) : std::abs(a));
133  return std::abs(a - b) <= smaller_magnitude * epsilon;
134  }
135 
136  // conversion from the ref element [0,1] (used by ESEAS) to the ref element [-1,1] (used by Intrepid2)
137  KOKKOS_INLINE_FUNCTION double fromZeroOne(double x_zero_one)
138  {
139  return x_zero_one * 2.0 - 1.0;
140  }
141 
142  // conversion from the ref element [-1,1] (used by Intrepid2) to the ref element [0,1] (used by ESEAS)
143  KOKKOS_INLINE_FUNCTION double toZeroOne(double x_minus_one_one)
144  {
145  return (x_minus_one_one + 1.0) / 2.0;
146  }
147 
148  // conversion from the ref element [0,1] (used by ESEAS) to the ref element [-1,1] (used by Intrepid2)
149  KOKKOS_INLINE_FUNCTION double fromZeroOne_dx(double dx_zero_one)
150  {
151  return dx_zero_one / 2.0;
152  }
153 
154  // conversion from the ref element [-1,1] (used by Intrepid2) to the ref element [0,1] (used by ESEAS)
155  KOKKOS_INLINE_FUNCTION double toZeroOne_dx(double dx_minus_one_one)
156  {
157  return dx_minus_one_one * 2.0;
158  }
159 
160  template<class DeviceViewType>
161  typename DeviceViewType::HostMirror getHostCopy(const DeviceViewType &deviceView)
162  {
163  typename DeviceViewType::HostMirror hostView = Kokkos::create_mirror(deviceView);
164  Kokkos::deep_copy(hostView, deviceView);
165  return hostView;
166  }
167 
168  template<class BasisFamily>
169  inline Teuchos::RCP< Intrepid2::Basis<DefaultTestDeviceType,double,double> > getBasisUsingFamily(shards::CellTopology cellTopo, Intrepid2::EFunctionSpace fs,
170  int polyOrder_x, int polyOrder_y=-1, int polyOrder_z = -1)
171  {
172  using BasisPtr = typename BasisFamily::BasisPtr;
173 
174  BasisPtr basis;
175  using namespace Intrepid2;
176 
177  if (cellTopo.getBaseKey() == shards::Line<>::key)
178  {
179  basis = getLineBasis<BasisFamily>(fs,polyOrder_x);
180  }
181  else if (cellTopo.getBaseKey() == shards::Quadrilateral<>::key)
182  {
183  INTREPID2_TEST_FOR_EXCEPTION(polyOrder_y < 0, std::invalid_argument, "polyOrder_y must be specified");
184  basis = getQuadrilateralBasis<BasisFamily>(fs,polyOrder_x,polyOrder_y);
185  }
186  else if (cellTopo.getBaseKey() == shards::Triangle<>::key)
187  {
188  basis = getTriangleBasis<BasisFamily>(fs,polyOrder_x);
189  }
190  else if (cellTopo.getBaseKey() == shards::Hexahedron<>::key)
191  {
192  INTREPID2_TEST_FOR_EXCEPTION(polyOrder_y < 0, std::invalid_argument, "polyOrder_y must be specified");
193  INTREPID2_TEST_FOR_EXCEPTION(polyOrder_z < 0, std::invalid_argument, "polyOrder_z must be specified");
194  basis = getHexahedronBasis<BasisFamily>(fs,polyOrder_x,polyOrder_y,polyOrder_z);
195  }
196  else if (cellTopo.getBaseKey() == shards::Tetrahedron<>::key)
197  {
198  basis = getTetrahedronBasis<BasisFamily>(fs, polyOrder_x);
199  }
200  else
201  {
202  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported cell topology");
203  }
204  return basis;
205  }
206 
207  template<bool defineVertexFunctions>
208  inline Teuchos::RCP< Intrepid2::Basis<DefaultTestDeviceType,double,double> > getHierarchicalBasis(shards::CellTopology cellTopo, Intrepid2::EFunctionSpace fs,
209  int polyOrder_x, int polyOrder_y=-1, int polyOrder_z = -1)
210  {
211  using DeviceType = DefaultTestDeviceType;
212  using Scalar = double;
213  using namespace Intrepid2;
214 
219 
221 
222  return getBasisUsingFamily<BasisFamily>(cellTopo, fs, polyOrder_x, polyOrder_y, polyOrder_z);
223  }
224 
225  template<typename ValueType, typename DeviceType, class ... DimArgs>
226  inline ViewType<ValueType,DeviceType> getView(const std::string &label, DimArgs... dims)
227  {
228  const bool allocateFadStorage = !std::is_pod<ValueType>::value;
229  if (!allocateFadStorage)
230  {
231  return ViewType<ValueType,DeviceType>(label,dims...);
232  }
233  else
234  {
235  return ViewType<ValueType,DeviceType>(label,dims...,MAX_FAD_DERIVATIVES_FOR_TESTS+1);
236  }
237  }
238 
239  // this method is to allow us to switch tests over incrementally; should collapse with ViewType once everything has been switched
240  template<typename ValueType, class ... DimArgs>
241  inline FixedRankViewType< typename RankExpander<ValueType, sizeof...(DimArgs) >::value_type, DefaultTestDeviceType > getFixedRankView(const std::string &label, DimArgs... dims)
242  {
243  const bool allocateFadStorage = !std::is_pod<ValueType>::value;
244  using value_type = typename RankExpander<ValueType, sizeof...(dims) >::value_type;
245  if (!allocateFadStorage)
246  {
247  return FixedRankViewType<value_type,DefaultTestDeviceType>(label,dims...);
248  }
249  else
250  {
251  return FixedRankViewType<value_type,DefaultTestDeviceType>(label,dims...,MAX_FAD_DERIVATIVES_FOR_TESTS+1);
252  }
253  }
254 
261  template <typename PointValueType, typename DeviceType>
262  inline ViewType<PointValueType,DeviceType> getInputPointsView(shards::CellTopology &cellTopo, int numPoints_1D)
263  {
264  const ordinal_type order = numPoints_1D - 1;
265  ordinal_type numPoints = PointTools::getLatticeSize(cellTopo, order);
266  ordinal_type spaceDim = cellTopo.getDimension();
267 
268  ViewType<PointValueType,DeviceType> inputPoints = getView<PointValueType,DeviceType>("input points",numPoints,spaceDim);
269  PointTools::getLattice(inputPoints, cellTopo, order, 0, POINTTYPE_EQUISPACED );
270 
271  return inputPoints;
272  }
273 
274  template<typename OutputValueType, typename DeviceType>
275  inline ViewType<OutputValueType,DeviceType> getOutputView(Intrepid2::EFunctionSpace fs, Intrepid2::EOperator op, int basisCardinality, int numPoints, int spaceDim)
276  {
277  switch (fs) {
278  case Intrepid2::FUNCTION_SPACE_HGRAD:
279  switch (op) {
280  case Intrepid2::OPERATOR_VALUE:
281  return getView<OutputValueType,DeviceType>("H^1 value output",basisCardinality,numPoints);
282  case Intrepid2::OPERATOR_GRAD:
283  return getView<OutputValueType,DeviceType>("H^1 derivative output",basisCardinality,numPoints,spaceDim);
284  case Intrepid2::OPERATOR_D1:
285  case Intrepid2::OPERATOR_D2:
286  case Intrepid2::OPERATOR_D3:
287  case Intrepid2::OPERATOR_D4:
288  case Intrepid2::OPERATOR_D5:
289  case Intrepid2::OPERATOR_D6:
290  case Intrepid2::OPERATOR_D7:
291  case Intrepid2::OPERATOR_D8:
292  case Intrepid2::OPERATOR_D9:
293  case Intrepid2::OPERATOR_D10:
294  {
295  const auto dkcard = getDkCardinality(op, spaceDim);
296  return getView<OutputValueType,DeviceType>("H^1 derivative output",basisCardinality,numPoints,dkcard);
297  }
298  default:
299  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported op/fs combination");
300  }
301  case Intrepid2::FUNCTION_SPACE_HCURL:
302  switch (op) {
303  case Intrepid2::OPERATOR_VALUE:
304  return getView<OutputValueType,DeviceType>("H(curl) value output",basisCardinality,numPoints,spaceDim);
305  case Intrepid2::OPERATOR_CURL:
306  if (spaceDim == 2)
307  return getView<OutputValueType,DeviceType>("H(curl) derivative output",basisCardinality,numPoints);
308  else if (spaceDim == 3)
309  return getView<OutputValueType,DeviceType>("H(curl) derivative output",basisCardinality,numPoints,spaceDim);
310  default:
311  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported op/fs combination");
312  }
313  case Intrepid2::FUNCTION_SPACE_HDIV:
314  switch (op) {
315  case Intrepid2::OPERATOR_VALUE:
316  return getView<OutputValueType,DeviceType>("H(div) value output",basisCardinality,numPoints,spaceDim);
317  case Intrepid2::OPERATOR_DIV:
318  return getView<OutputValueType,DeviceType>("H(div) derivative output",basisCardinality,numPoints);
319  default:
320  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported op/fs combination");
321  }
322 
323  case Intrepid2::FUNCTION_SPACE_HVOL:
324  switch (op) {
325  case Intrepid2::OPERATOR_VALUE:
326  return getView<OutputValueType,DeviceType>("H(vol) value output",basisCardinality,numPoints);
327  case Intrepid2::OPERATOR_D1:
328  case Intrepid2::OPERATOR_D2:
329  case Intrepid2::OPERATOR_D3:
330  case Intrepid2::OPERATOR_D4:
331  case Intrepid2::OPERATOR_D5:
332  case Intrepid2::OPERATOR_D6:
333  case Intrepid2::OPERATOR_D7:
334  case Intrepid2::OPERATOR_D8:
335  case Intrepid2::OPERATOR_D9:
336  case Intrepid2::OPERATOR_D10:
337  {
338  const auto dkcard = getDkCardinality(op, spaceDim);
339  return getView<OutputValueType,DeviceType>("H(vol) derivative output",basisCardinality,numPoints,dkcard);
340  }
341  default:
342  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported op/fs combination");
343  }
344  default:
345  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported op/fs combination");
346  }
347  }
348 
349  // ! This returns a vector whose entries are vector<int>s containing 1-3 polynomial orders from 1 up to and including those specified
350  // ! Intended for testing bases that support anisotropic polynomial degree, such as the hierarchical bases
351  inline std::vector< std::vector<int> > getBasisTestCasesUpToDegree(int spaceDim, int minDegree, int polyOrder_x, int polyOrder_y=-1, int polyOrder_z = -1)
352  {
353  std::vector<int> degrees(spaceDim);
354  degrees[0] = polyOrder_x;
355  if (spaceDim > 1) degrees[1] = polyOrder_y;
356  if (spaceDim > 2) degrees[2] = polyOrder_z;
357 
358  int numCases = degrees[0];
359  for (unsigned d=1; d<degrees.size(); d++)
360  {
361  numCases = numCases * (degrees[d] + 1 - minDegree);
362  }
363  std::vector< std::vector<int> > subBasisDegreeTestCases(numCases);
364  for (int caseOrdinal=0; caseOrdinal<numCases; caseOrdinal++)
365  {
366  std::vector<int> subBasisDegrees(degrees.size());
367  int caseRemainder = caseOrdinal;
368  for (int d=degrees.size()-1; d>=0; d--)
369  {
370  int subBasisDegree = caseRemainder % (degrees[d] + 1 - minDegree);
371  caseRemainder = caseRemainder / (degrees[d] + 1 - minDegree);
372  subBasisDegrees[d] = subBasisDegree + minDegree;
373  }
374  subBasisDegreeTestCases[caseOrdinal] = subBasisDegrees;
375  }
376  return subBasisDegreeTestCases;
377  }
378 
380  template<class Functor, class Scalar, int rank>
381  typename ViewType<Scalar,DefaultTestDeviceType>::HostMirror copyFunctorToHostView(const Functor &deviceFunctor)
382  {
383  INTREPID2_TEST_FOR_EXCEPTION(rank != getFunctorRank(deviceFunctor), std::invalid_argument, "functor rank must match the template argument");
384 
385  using DeviceType = DefaultTestDeviceType;
386  ViewType<Scalar,DeviceType> view;
387  const std::string label = "functor copy";
388  const auto &f = deviceFunctor;
389  switch (rank)
390  {
391  case 0:
392  view = getView<Scalar,DeviceType>(label);
393  break;
394  case 1:
395  view = getView<Scalar,DeviceType>(label, f.extent_int(0));
396  break;
397  case 2:
398  view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1));
399  break;
400  case 3:
401  view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2));
402  break;
403  case 4:
404  view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2), f.extent_int(3));
405  break;
406  case 5:
407  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));
408  break;
409  case 6:
410  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));
411  break;
412  case 7:
413  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));
414  break;
415  default:
416  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported functor rank");
417  }
418 
419  int entryCount = view.size();
420 
421  using ExecutionSpace = typename ViewType<Scalar,DeviceType>::execution_space;
422 
423  using ViewIteratorScalar = Intrepid2::ViewIterator<ViewType<Scalar,DeviceType>, Scalar>;
424  using FunctorIteratorScalar = FunctorIterator<Functor, Scalar, rank>;
425 
426  Kokkos::RangePolicy < ExecutionSpace > policy(0,entryCount);
427  Kokkos::parallel_for( policy,
428  KOKKOS_LAMBDA (const int &enumerationIndex )
429  {
430  ViewIteratorScalar vi(view);
431  FunctorIteratorScalar fi(f);
432 
433  vi.setEnumerationIndex(enumerationIndex);
434  fi.setEnumerationIndex(enumerationIndex);
435 
436  vi.set(fi.get());
437  }
438  );
439 
440  auto hostView = Kokkos::create_mirror_view(view);
441  Kokkos::deep_copy(hostView, view);
442  return hostView;
443  }
444 
445  template<class FunctorType, typename Scalar, int rank>
446  void printFunctor(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
447  {
448  auto functorHostCopy = copyFunctorToHostView<FunctorType, Scalar, rank>(functor);
449 
450  std::string name = (functorName == "") ? "Functor" : functorName;
451 
452  out << "\n******** " << name << " (rank " << rank << ") ********\n";
453  out << "dimensions: (";
454  for (int r=0; r<rank; r++)
455  {
456  out << functor.extent_int(r);
457  if (r<rank-1) out << ",";
458  }
459  out << ")\n";
460 
462 
463  bool moreEntries = true;
464  while (moreEntries)
465  {
466  Scalar value = vi.get();
467 
468  auto location = vi.getLocation();
469  out << functorName << "(";
470  for (ordinal_type i=0; i<rank; i++)
471  {
472  out << location[i];
473  if (i<rank-1)
474  {
475  out << ",";
476  }
477  }
478  out << ") " << value << std::endl;
479 
480  moreEntries = (vi.increment() != -1);
481  }
482  out << "\n";
483  }
484 
485  template<class FunctorType>
486  void printFunctor1(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
487  {
488  using Scalar = typename std::remove_reference<decltype(functor(0))>::type;
489  printFunctor<FunctorType, Scalar, 1>(functor, out, functorName);
490  }
491 
492  template<class FunctorType>
493  void printFunctor2(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
494  {
495  using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0))>::type>::type;
496  printFunctor<FunctorType, Scalar, 2>(functor, out, functorName);
497  }
498 
499  template<class FunctorType>
500  void printFunctor3(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
501  {
502  using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0,0))>::type>::type;
503  printFunctor<FunctorType, Scalar, 3>(functor, out, functorName);
504  }
505 
506  template<class FunctorType>
507  void printFunctor4(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
508  {
509  using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0,0,0))>::type>::type;
510  printFunctor<FunctorType, Scalar, 4>(functor, out, functorName);
511  }
512 
513  template<class FunctorType>
514  void printFunctor5(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
515  {
516  using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0,0,0,0))>::type>::type;
517  printFunctor<FunctorType, Scalar, 5>(functor, out, functorName);
518  }
519 
520  template<class FunctorType>
521  void printFunctor6(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
522  {
523  using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0,0,0,0,0))>::type>::type;
524  printFunctor<FunctorType, Scalar, 6>(functor, out, functorName);
525  }
526 
527  template<class FunctorType>
528  void printFunctor7(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
529  {
530  using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0,0,0,0,0,0))>::type>::type;
531  printFunctor<FunctorType, Scalar, 7>(functor, out, functorName);
532  }
533 
534  template<class View>
535  void printView(const View &view, std::ostream &out, const std::string &viewName = "")
536  {
537  using Scalar = typename View::value_type;
538  using HostView = typename View::HostMirror;
539  using HostViewIteratorScalar = Intrepid2::ViewIterator<HostView, Scalar>;
540 
541  auto hostView = getHostCopy(view);
542 
543  HostViewIteratorScalar vi(hostView);
544 
545  bool moreEntries = (vi.nextIncrementRank() != -1);
546  while (moreEntries)
547  {
548  Scalar value = vi.get();
549 
550  auto location = vi.getLocation();
551  out << viewName << "(";
552  for (unsigned i=0; i<getFunctorRank(view); i++)
553  {
554  out << location[i];
555  if (i<getFunctorRank(view)-1)
556  {
557  out << ",";
558  }
559  }
560  out << ") " << value << std::endl;
561 
562  moreEntries = (vi.increment() != -1);
563  }
564  }
565 
566  template <class FunctorType1, class FunctorType2, int rank, typename Scalar=typename FunctorType1::value_type, class ExecutionSpace = typename FunctorType1::execution_space>
567  typename std::enable_if< !supports_rank<FunctorType1,rank>::value, void >::type
568  testFloatingEquality(const FunctorType1 &functor1, const FunctorType2 &functor2, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
569  std::string functor1Name = "Functor 1", std::string functor2Name = "Functor 2")
570  {
571  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "testFloatingEquality() called on FunctorType1 that does not support the specified rank.\n");
572  }
573 
575  template <class FunctorType1, class FunctorType2, int rank, typename Scalar=typename FunctorType1::value_type, class ExecutionSpace = typename FunctorType1::execution_space>
576  typename std::enable_if< supports_rank<FunctorType1,rank>::value, void >::type
577  testFloatingEquality(const FunctorType1 &functor1, const FunctorType2 &functor2, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
578  std::string functor1Name = "Functor 1", std::string functor2Name = "Functor 2")
579  {
580  static_assert( supports_rank<FunctorType2,rank>::value, "Both Functor1 and Functor2 must support the specified rank through operator().");
581  using Functor1IteratorScalar = FunctorIterator<FunctorType1, Scalar, rank>;
582  using Functor2IteratorScalar = FunctorIterator<FunctorType2, Scalar, rank>;
583 
584  // check that rank/size match
585  TEUCHOS_TEST_FOR_EXCEPTION(getFunctorRank(functor1) != rank, std::invalid_argument, "functor1 and functor2 must agree in rank"); // Kokkos::View does not provide rank() method; getFunctorRank() provides common interface
586  TEUCHOS_TEST_FOR_EXCEPTION(getFunctorRank(functor2) != rank, std::invalid_argument, "functor1 and functor2 must agree in rank"); // Kokkos::View does not provide rank() method; getFunctorRank() provides common interface
587 
588  int entryCount = 1;
589  for (unsigned i=0; i<getFunctorRank(functor1); i++)
590  {
591  TEUCHOS_TEST_FOR_EXCEPTION(functor1.extent_int(i) != functor2.extent_int(i), std::invalid_argument,
592  "functor1 and functor2 must agree in size in each dimension; functor1 has extent "
593  + std::to_string(functor1.extent_int(i)) + " in dimension " + std::to_string(i)
594  + "; functor2 has extent " + std::to_string(functor2.extent_int(i)) );
595  entryCount *= functor1.extent_int(i);
596  }
597  if (entryCount == 0) return; // nothing to test
598 
599  ViewType<bool,ExecutionSpace> valuesMatch = getView<bool,ExecutionSpace>("valuesMatch", entryCount);
600 
601  Kokkos::RangePolicy < ExecutionSpace > policy(0,entryCount);
602  Kokkos::parallel_for( policy,
603  KOKKOS_LAMBDA (const int &enumerationIndex )
604  {
605  Functor1IteratorScalar vi1(functor1);
606  Functor2IteratorScalar vi2(functor2);
607 
608  vi1.setEnumerationIndex(enumerationIndex);
609  vi2.setEnumerationIndex(enumerationIndex);
610 
611  const Scalar & value1 = vi1.get();
612  const Scalar & value2 = vi2.get();
613 
614  bool errMeetsTol = errMeetsAbsAndRelTol(value1, value2, relTol, absTol);
615  valuesMatch(enumerationIndex) = errMeetsTol;
616  }
617  );
618 
619  int failureCount = 0;
620  Kokkos::RangePolicy<ExecutionSpace > reducePolicy(0, entryCount);
621  Kokkos::parallel_reduce( reducePolicy,
622  KOKKOS_LAMBDA( const int &enumerationIndex, int &reducedValue )
623  {
624  reducedValue += valuesMatch(enumerationIndex) ? 0 : 1;
625  }, failureCount);
626 
627  const bool allValuesMatch = (failureCount == 0);
628  success = success && allValuesMatch;
629 
630  if (!allValuesMatch)
631  {
632  // copy functors to host views
633  auto functor1CopyHost = copyFunctorToHostView<FunctorType1,Scalar,rank>(functor1);
634  auto functor2CopyHost = copyFunctorToHostView<FunctorType2,Scalar,rank>(functor2);
635 
636  auto valuesMatchHost = getHostCopy(valuesMatch);
637 
641 
642  bool moreEntries = true;
643  while (moreEntries)
644  {
645  bool errMeetsTol = viMatch.get();
646 
647  if (!errMeetsTol)
648  {
649  const Scalar value1 = vi1.get();
650  const Scalar value2 = vi2.get();
651 
652  success = false;
653  auto location = vi1.getLocation();
654  out << "At location (";
655  for (unsigned i=0; i<getFunctorRank(functor1); i++)
656  {
657  out << location[i];
658  if (i<getFunctorRank(functor1)-1)
659  {
660  out << ",";
661  }
662  }
663  out << ") " << functor1Name << " value != " << functor2Name << " value (";
664  out << value1 << " != " << value2 << "); diff is " << std::abs(value1-value2) << std::endl;
665  }
666 
667  moreEntries = (vi1.increment() != -1);
668  moreEntries = moreEntries && (vi2.increment() != -1);
669  moreEntries = moreEntries && (viMatch.increment() != -1);
670  }
671  }
672  }
673 
674  template <class ViewType, class FunctorType>
675  void testFloatingEquality1(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
676  std::string view1Name = "View", std::string view2Name = "Functor")
677  {
678  testFloatingEquality<ViewType, FunctorType, 1>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
679  }
680 
681  template <class ViewType, class FunctorType>
682  void testFloatingEquality2(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
683  std::string view1Name = "View", std::string view2Name = "Functor")
684  {
685  testFloatingEquality<ViewType, FunctorType, 2>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
686  }
687 
688  template <class ViewType, class FunctorType>
689  void testFloatingEquality3(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
690  std::string view1Name = "View", std::string view2Name = "Functor")
691  {
692  testFloatingEquality<ViewType, FunctorType, 3>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
693  }
694 
695  template <class ViewType, class FunctorType>
696  void testFloatingEquality4(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
697  std::string view1Name = "View", std::string view2Name = "Functor")
698  {
699  testFloatingEquality<ViewType, FunctorType, 4>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
700  }
701 
702  template <class ViewType, class FunctorType>
703  void testFloatingEquality5(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
704  std::string view1Name = "View", std::string view2Name = "Functor")
705  {
706  testFloatingEquality<ViewType, FunctorType, 5>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
707  }
708 
709  template <class ViewType, class FunctorType>
710  void testFloatingEquality6(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
711  std::string view1Name = "View", std::string view2Name = "Functor")
712  {
713  testFloatingEquality<ViewType, FunctorType, 6>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
714  }
715 
716  template <class ViewType, class FunctorType>
717  void testFloatingEquality7(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
718  std::string view1Name = "View", std::string view2Name = "Functor")
719  {
720  testFloatingEquality<ViewType, FunctorType, 7>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
721  }
722 
723  template <class ViewType1, class ViewType2>
724  void testViewFloatingEquality(const ViewType1 &view1, const ViewType2 &view2, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
725  std::string view1Name = "View 1", std::string view2Name = "View 2")
726  {
727  // check that rank/size match
728  TEUCHOS_TEST_FOR_EXCEPTION(view1.rank() != view2.rank(), std::invalid_argument, "views must agree in rank");
729  for (unsigned i=0; i<view1.rank(); i++)
730  {
731  TEUCHOS_TEST_FOR_EXCEPTION(view1.extent_int(i) != view2.extent_int(i), std::invalid_argument, "views must agree in size in each dimension");
732  }
733 
734  if (view1.size() == 0) return; // nothing to test
735 
736  const int rank = view1.rank();
737  switch (rank)
738  {
739  case 0: testFloatingEquality<ViewType1, ViewType2, 0>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
740  case 1: testFloatingEquality<ViewType1, ViewType2, 1>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
741  case 2: testFloatingEquality<ViewType1, ViewType2, 2>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
742  case 3: testFloatingEquality<ViewType1, ViewType2, 3>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
743  case 4: testFloatingEquality<ViewType1, ViewType2, 4>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
744  case 5: testFloatingEquality<ViewType1, ViewType2, 5>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
745  case 6: testFloatingEquality<ViewType1, ViewType2, 6>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
746  case 7: testFloatingEquality<ViewType1, ViewType2, 7>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
747  default: INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported rank");
748  }
749  }
750 
751 } // namespace Intrepid2
752 
753 #ifdef HAVE_INTREPID2_SACADO
754 using Sacado_Fad_DFadType = Sacado::Fad::DFad<double>;
755 #define INTREPID2_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
756 \
757 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( GROUP_NAME, TEST_NAME, double ) \
758 \
759 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( GROUP_NAME, TEST_NAME, Sacado_Fad_DFadType ) \
760 
761 #define INTREPID2_OUTPUTSCALAR_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
762 \
763 TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( GROUP_NAME, TEST_NAME, double, double ) \
764 \
765 TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( GROUP_NAME, TEST_NAME, Sacado_Fad_DFadType, Sacado_Fad_DFadType ) \
766 
767 #else
768 #define INTREPID2_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
769 \
770 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( GROUP_NAME, TEST_NAME, double ) \
771 
772 #define INTREPID2_OUTPUTSCALAR_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
773 \
774 TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( GROUP_NAME, TEST_NAME, double, double ) \
775 
776 #endif
777 
778 #endif /* Intrepid2_TestUtils_h */
enable_if_t< M==0 > KOKKOS_INLINE_FUNCTION get() const
Stateless classes that act as factories for two families of hierarchical bases. HierarchicalBasisFami...
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
KOKKOS_INLINE_FUNCTION int increment()
Stateless class representing a family of basis functions, templated on H(vol) and H(grad) on the line...
KOKKOS_INLINE_FUNCTION Kokkos::Array< int, 7 > & getLocation()
KOKKOS_INLINE_FUNCTION ScalarType get()
SFINAE helper to detect whether a type supports a rank-integral-argument operator().
Header function for Intrepid2::Util class and other utility functions.
ViewType< Scalar, DefaultTestDeviceType >::HostMirror copyFunctorToHostView(const Functor &deviceFunctor)
Copy the values for the specified functor.
Helper to get Scalar[*+] where the number of *&#39;s matches the given rank.
KOKKOS_INLINE_FUNCTION ordinal_type getDkCardinality(const EOperator operatorType, const ordinal_type spaceDim)
Returns multiplicities of dx, dy, and dz based on the enumeration of the partial derivative, its order and the space dimension. Inverse of the getDkEnumeration() method.
Basis defining Legendre basis on the line, a polynomial subspace of L^2 (a.k.a. H(vol)) on the line...
constexpr int MAX_FAD_DERIVATIVES_FOR_TESTS
Maximum number of derivatives to track for Fad types in tests.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line...
ViewType< PointValueType, DeviceType > getInputPointsView(shards::CellTopology &cellTopo, int numPoints_1D)
Returns a DynRankView containing regularly-spaced points on the specified cell topology.
essentially, a read-only variant of ViewIterator, for a general functor (extent_int() and rank() supp...
Defines the Intrepid2::FunctorIterator class, as well as the Intrepid2::functor_returns_ref SFINAE he...
A family of basis functions, constructed from H(vol) and H(grad) bases on the line.
KOKKOS_INLINE_FUNCTION int increment()
Header file to include all Sacado headers that are required if using Intrepid2 with Sacado types...
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.
const Teuchos::ScalarTraits< Scalar >::magnitudeType smallNumber()
Use Teuchos small number determination on host; pass this to Intrepid2::relErr() on device...
static void getLattice(Kokkos::DynRankView< pointValueType, pointProperties... > points, const shards::CellTopology cellType, const ordinal_type order, const ordinal_type offset=0, const EPointType pointType=POINTTYPE_EQUISPACED)
Computes a lattice of points of a given order on a reference simplex, quadrilateral or hexahedron (cu...
A helper class that allows iteration over some part of a Kokkos View, while allowing the calling code...
enable_if_t< has_rank_method< Functor >::value, unsigned > KOKKOS_INLINE_FUNCTION getFunctorRank(const Functor &functor)
static ordinal_type getLatticeSize(const shards::CellTopology cellType, const ordinal_type order, const ordinal_type offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...
typename Kokkos::DefaultExecutionSpace::device_type DefaultTestDeviceType
Default Kokkos::Device to use for tests; depends on platform.
Header file for the abstract base class Intrepid2::Basis.
Header file for Intrepid2::PointTools class to provide utilities for barycentric coordinates, equispaced lattices, and warp-blend point distrubtions.