Quadmath

src/amdis/DataTransfer.inc.hpp

@@ -365,7 +365,7 @@ restrictLocal(Element const& father, NodeElementData& fatherDOFs, Trafo const& t
         for (std::size_t i = 0; i < shapeValues.size(); ++i)
           y += shapeValues[i] * childDOFs[i];
 
-        fatherDOFsTemp_[currentDOF] = y;
+        fatherDOFsTemp_[currentDOF] = T(y);
         finishedDOFs_[currentDOF++] = true;
         return y;
       }

src/amdis/common/CMakeLists.txt

@@ -27,6 +27,7 @@ install(FILES
     MultiTypeVector.hpp
     Range.hpp
     Resize.hpp
+    QuadMath.hpp
     StaticSize.hpp
     String.hpp
     SwitchCases.hpp

src/amdis/common/FieldMatVec.inc.hpp

@@ -346,13 +346,15 @@ auto two_norm(T const& x)
 template <class T, int N>
 auto two_norm(FieldVector<T, N> const& x)
 {
-  return std::sqrt(unary_dot(x));
+  using std::sqrt;
+  return sqrt(unary_dot(x));
 }
 
 template <class T, int N>
 auto two_norm(FieldMatrix<T, 1, N> const& x)
 {
-  return std::sqrt(unary_dot(x));
+  using std::sqrt;
+  return sqrt(unary_dot(x));
 }
 
 /** \ingroup vector_norms

src/amdis/common/QuadMath.hpp

+#pragma once
+
+#if HAVE_QUADMATH
+#include <functional>
+#include <type_traits>
+
+#include <dune/common/hash.hh>
+#include <dune/common/quadmath.hh>
+
+#include <amdis/common/DerivativeTraits.hpp>
+#include <amdis/common/StaticSize.hpp>
+#include <amdis/common/ValueCategory.hpp>
+#include <amdis/gridfunctions/ConstantGridFunction.hpp>
+
+namespace std
+{
+  template <class T>
+  struct common_type<Dune::Float128, T>
+  {
+    using type = Dune::Float128;
+  };
+
+  template <class T>
+  struct common_type<T, Dune::Float128>
+  {
+    using type = Dune::Float128;
+  };
+
+  template <>
+  struct common_type<Dune::Float128, Dune::Float128>
+  {
+    using type = Dune::Float128;
+  };
+
+  template <>
+  struct hash<Dune::Float128>
+  {
+    typedef Dune::Float128 argument_type;
+    typedef std::size_t result_type;
+
+    std::size_t operator()(const Dune::Float128& arg) const
+    {
+      hash<long double> hasher_ld;
+      return hasher_ld((long double)(arg));
+    }
+  };
+
+  template <>
+  struct hash<const Dune::Float128>
+      : public hash<Dune::Float128>
+  {};
+
+} // end namespace std
+
+
+namespace Dune
+{
+  namespace Impl
+  {
+    // specialization for float arguments due to ambiguity
+    template <class T,
+      std::enable_if_t<not std::is_integral<T>::value && std::is_arithmetic<T>::value, int> = 0>
+    inline Float128 pow(const Float128& x, const T& p)
+    {
+      return powq(float128_t(x), float128_t(p));
+    }
+
+  } // end namespace Impl
+} // end namespace Dune
+
+
+namespace AMDiS
+{
+  namespace Concepts
+  {
+    namespace Definition
+    {
+      template <>
+      struct ConstantToGridFunction<Dune::Float128>
+        : std::true_type {};
+
+    } // end namespace Definition
+  } // end namespace Concepts
+
+  namespace Impl
+  {
+    template <>
+    struct SizeImpl<Dune::Float128>
+        : std::integral_constant<std::size_t, 1> {};
+
+    template <>
+    struct RowsImpl<Dune::Float128>
+        : std::integral_constant<std::size_t, 1> {};
+
+    template <>
+    struct ColsImpl<Dune::Float128>
+        : std::integral_constant<std::size_t, 1> {};
+
+  } // end namespace Impl
+
+  template <class K, int N>
+  struct DerivativeTraits<Dune::Float128(Dune::FieldVector<K,N>), tag::gradient>
+  {
+    using Range = Dune::FieldVector<Dune::Float128,N>;
+  };
+
+  template <>
+  struct ValueCategory<Dune::Float128>
+  {
+    using type = tag::scalar;
+  };
+
+} // end namespace AMDiS
+
+#endif // HAVE_QUADMATH

src/amdis/linearalgebra/petsc/Traits.hpp

@@ -16,7 +16,7 @@ namespace AMDiS
   template <class Basis, class T = double>
   struct BackendTraits
   {
-    static_assert(std::is_convertible<T,PetscScalar>::value, "");
+    static_assert(std::is_same<T,PetscScalar>::value, "");
 
     using Mat = ::Mat;
     using Vec = ::Vec;

test/ProblemStatTest.cpp

@@ -2,20 +2,33 @@
 #include <amdis/AdaptStationary.hpp>
 #include <amdis/LocalOperators.hpp>
 #include <amdis/ProblemStat.hpp>
+#include <amdis/common/QuadMath.hpp>
 
-using namespace AMDiS;
-
-using Grid = Dune::YaspGrid<2>;
+#include <dune/common/version.hh>
+#include <dune/grid/yaspgrid.hh>
 
-template <class T>
-struct Param
-    : DefaultBasisCreator<Grid, Impl::LagrangePreBasisCreator<1,1>, T>
-{};
+using namespace AMDiS;
 
 template <class T>
 void test()
 {
-  ProblemStat<Param<T>> prob("ellipt");
+  // use T as coordinate type
+  using Grid = Dune::YaspGrid<2, Dune::EquidistantCoordinates<T,2>>;
+  Grid grid({T(1), T(1)}, {8,8});
+
+  // use T as range type for basis
+  using namespace Dune::Functions::BasisFactory;
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+  auto basis = makeBasis(grid.leafGridView(), power<1>(lagrange<1>(), flatLexicographic()));
+#else
+  auto basis = makeBasis(grid.leafGridView(), power<1>(lagrange<1,T>(), flatLexicographic()));
+#endif
+  using Basis = decltype(basis);
+
+  // use T as coefficient type
+  using Param = DefaultProblemTraits<Basis, T>;
+
+  ProblemStat<Param> prob("ellipt", grid, basis);
   prob.initialize(INIT_ALL);
   prob.boundaryManager()->setBoxBoundary({1,1,2,2});
 
@@ -45,5 +58,9 @@ int main(int argc, char** argv)
   test<double>();
   test<long double>();
 
+#if HAVE_QUADMATH
+  test<Dune::Float128>();
+#endif
+
   return 0;
 }