MultiType Vector/Matrix and some tests added

CMakeLists.txt

@@ -21,6 +21,7 @@ dune_enable_all_packages(MODULE_LIBRARIES amdis)
 
 include(AmdisMacros)
 add_subdirectory("src")
+add_subdirectory("test")
 add_subdirectory("dune")
 add_subdirectory("doc")
 add_subdirectory("cmake/modules")

dune/amdis/FileWriter.hpp

@@ -91,6 +91,8 @@ namespace AMDiS
     /// Implements \ref FileWriterInterface::writeFiles
     virtual void writeFiles(AdaptInfo& adaptInfo, bool force) override
     {
+      if (!filesystem::exists(dir_))
+        error_exit("Output directory '",dir_,"' does not exist!");
       if (vtkSeqWriter_)
         vtkSeqWriter_->write(adaptInfo.getTime(), mode_);
     }

dune/amdis/Mesh.hpp

@@ -146,7 +146,7 @@ namespace AMDiS
 
       msg("L = ", L);
 
-      auto s = Dune::filledArray<std::size_t(dim)>(2); // number of cells on coarse mesh in each direction
+      auto s = Dune::filledArray<std::size_t(dim)>(1); // number of cells on coarse mesh in each direction
       Parameters::get(meshName + "->num cells", s);
 
       // TODO: add more parameters for yasp-grid (see constructor)
@@ -168,7 +168,7 @@ namespace AMDiS
       Parameters::get(meshName + "->min corner", lowerleft);
       Parameters::get(meshName + "->max corner", upperright);
 
-      auto s = Dune::filledArray<std::size_t(dim)>(2); // number of cells on coarse mesh in each direction
+      auto s = Dune::filledArray<std::size_t(dim)>(1); // number of cells on coarse mesh in each direction
       Parameters::get(meshName + "->num cells", s);
 
       // TODO: add more parameters for yasp-grid (see constructor)

dune/amdis/assembler/FirstOrderTestPartialTrial.hpp

@@ -59,6 +59,7 @@ namespace AMDiS
 
         // The transposed inverse Jacobian of the map from the reference element to the element
         const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+        assert(jacobian.M() == Context::dim);
 
         // The multiplicative factor in the integral transformation formula
         double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();

dune/amdis/assembler/ZeroOrderTestvecTrialvec.hpp

@@ -41,8 +41,8 @@ namespace AMDiS
                                 ElementMatrix& elementMatrix,
                                 RowNode const& rowNode,
                                 ColNode const& colNode,
-                                std::integral_constant<bool, sameFE>,
-                                std::integral_constant<bool, sameNode>)
+                                bool_t<sameFE>,
+                                bool_t<sameNode>)
     {
       static_assert(RowNode::isPower && ColNode::isPower,
         "Operator can be applied to Power-Nodes only.");
@@ -50,12 +50,10 @@ namespace AMDiS
       static_assert( (RowNode::CHILDREN == ColNode::CHILDREN), "" );
       // theoretically the restriction of the same nr of childs would not be necessary
 
-      using Category = ValueCategory_t<typename GridFct::Range>;
-      static_assert( std::is_same<Category, tag::scalar>::value || std::is_same<Category, tag::matrix>::value, "" );
-
+      using Category = ValueCategory_t<expr_value_type>;
       calcElementMatrix(context, quad, elementMatrix, rowNode, colNode,
-        std::integral_constant<bool, sameFE>{},
-        std::integral_constant<bool, sameNode>{},
+        bool_<sameFE>,
+        bool_<sameNode>,
         Category{});
     }
 
@@ -68,7 +66,7 @@ namespace AMDiS
                            ElementMatrix& elementMatrix,
                            RowNode const& rowNode,
                            ColNode const& colNode,
-                           std::integral_constant<bool, sameFE>,
+                           bool_t<sameFE>,
                            std::false_type /*sameNode*/,
                            tag::scalar)
     {
@@ -165,12 +163,12 @@ namespace AMDiS
                            ElementMatrix& elementMatrix,
                            RowNode const& rowNode,
                            ColNode const& colNode,
-                           std::integral_constant<bool, sameFE>,
-                           std::integral_constant<bool, sameNode>,
+                           bool_t<sameFE>,
+                           bool_t<sameNode>,
                            tag::matrix)
     {
       static const std::size_t CHILDREN = RowNode::CHILDREN;
-      static_assert( std::is_same<typename GridFct::Range, Dune::FieldMatrix<double, CHILDREN, CHILDREN>>::value, "" );
+      static_assert( std::is_same<expr_value_type, Dune::FieldMatrix<double, CHILDREN, CHILDREN>>::value, "" );
 
       auto const& rowLocalFE = rowNode.child(0).finiteElement();
       auto const& colLocalFE = colNode.child(0).finiteElement();
@@ -191,10 +189,10 @@ namespace AMDiS
           colLocalFE.localBasis().evaluateFunction(local, colShapeValues_);
 
         for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
-          double value = factor * (*rowShapeValues)[i];
+          double value0 = factor * (*rowShapeValues)[i];
 
           for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
-            value *= (*colShapeValues)[j];
+            double value = value0 * (*colShapeValues)[j];
             auto mat = exprValue * value;
 
             for (std::size_t k0 = 0; k0 < CHILDREN; ++k0) {

dune/amdis/common/CMakeLists.txt

@@ -11,7 +11,7 @@ install(FILES
     Loops.hpp
     Math.hpp
     Mpl.hpp
-    MultiTypeFieldVector.hpp
+    MultiTypeVector.hpp
     ScalarTypes.hpp
     Size.hpp
     TupleUtility.hpp

dune/amdis/common/Concepts.hpp

@@ -108,6 +108,19 @@ namespace AMDiS
           >;
       };
 
+      // idx[0]
+      struct MultiIndex
+      {
+        template <class MI>
+        auto requires_(MI&& idx) -> decltype(
+          Concepts::valid_expr(
+            idx[0],
+            idx.size(),
+            idx.max_size()
+            /* ,idx.resize() */
+          ));
+      };
+
     } // end namespace Definition
 #endif // DOXYGEN
 
@@ -166,6 +179,9 @@ namespace AMDiS
     template <class F, class... Args>
     constexpr bool Predicate = Functor<F, bool(Args...)>;
 
+    template <class MI>
+    constexpr bool MultiIndex = models<Definition::MultiIndex(MI)>;
+
     /** @} **/
 
   } // end namespace Concepts

dune/amdis/common/FieldMatVec.hpp

@@ -2,6 +2,7 @@
 
 #include <algorithm>
 
+#include <dune/common/diagonalmatrix.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 
@@ -371,20 +372,23 @@ namespace AMDiS
   }
 
 
-  template <class T, int M, int N>
-  FieldMatrix<T,M,N> operator*(T scalar, FieldMatrix<T, M, N> A)
+  template <class T, int M, int N, class S,
+    REQUIRES(Concepts::Arithmetic<S>) >
+  FieldMatrix<T,M,N> operator*(S scalar, FieldMatrix<T, M, N> A)
   {
     return A *= scalar;
   }
 
-  template <class T, int M, int N>
-  FieldMatrix<T,M,N> operator*(FieldMatrix<T, M, N> A, T scalar)
+  template <class T, int M, int N, class S,
+    REQUIRES(Concepts::Arithmetic<S>) >
+  FieldMatrix<T,M,N> operator*(FieldMatrix<T, M, N> A, S scalar)
   {
     return A *= scalar;
   }
 
-  template <class T, int M, int N>
-  FieldMatrix<T,M,N> operator/(FieldMatrix<T, M, N> A, T scalar)
+  template <class T, int M, int N, class S,
+    REQUIRES(Concepts::Arithmetic<S>) >
+  FieldMatrix<T,M,N> operator/(FieldMatrix<T, M, N> A, S scalar)
   {
     return A /= scalar;
   }

dune/amdis/common/FieldTraits.hpp

+#pragma once
+
+#include <dune/common/ftraits.hh>
+
+#include <dune/amdis/common/Utility.hpp>
+
+namespace AMDiS
+{
+  template <class... T>
+  struct CommonFieldTraits
+  {
+    using field_type = CommonType_t<typename Dune::FieldTraits<T>::field_type...>;
+    using real_type = CommonType_t<typename Dune::FieldTraits<T>::real_type...>;
+  };
+
+} // end namespace AMDiS

dune/amdis/common/Mpl.hpp

 #pragma once
 
 // std c++ headers
+#include <tuple>
 #include <type_traits>
 #include <utility>
 
@@ -152,8 +153,8 @@ namespace AMDiS
   constexpr bool_t<!B> operator!(bool_t<B>) { return {}; }
 
 
-  template <class T, T A, T B>
-  using is_equal = std::conditional_t<A == B, std::true_type, std::false_type>;
+  template <class T, T value0, T... values>
+  using is_equal = all_of_t<(value0 == values)...>;
 
   template <class T, class... Ts>
   using is_one_of = or_t<std::is_same<T, Ts>::value...>;

dune/amdis/common/MultiTypeMatrix.hpp

+#pragma once
+
+#include <tuple>
+
+#include <dune/amdis/common/Concepts.hpp>
+#include <dune/amdis/common/FieldTraits.hpp>
+#include <dune/amdis/common/Loops.hpp>
+#include <dune/amdis/common/Mpl.hpp>
+#include <dune/amdis/common/MultiTypeVector.hpp>
+#include <dune/amdis/common/Size.hpp>
+#include <dune/amdis/utility/MultiIndex.hpp>
+
+namespace AMDiS
+{
+  // forward declaration
+  template <class... Rows>
+  class MultiTypeMatrix;
+}
+
+namespace Dune
+{
+  template <class... Rows>
+  struct FieldTraits<AMDiS::MultiTypeMatrix<Rows...>>
+  {
+    using field_type = typename AMDiS::CommonFieldTraits<Rows...>::field_type;
+    using real_type = typename AMDiS::CommonFieldTraits<Rows...>::real_type;
+  };
+}
+
+
+namespace AMDiS
+{
+  // Rows should be of type MultiTypeVector
+  template <class... Rows>
+  class MultiTypeMatrix
+      : public std::tuple<Rows...>
+  {
+    using Self = MultiTypeMatrix;
+    using Super = std::tuple<Rows...>;
+
+    static_assert(is_equal<int, Rows::dimension...>::value,
+      "All columns must have the same length.");
+
+  public:
+    using field_type = typename Dune::FieldTraits<Self>::field_type;
+    using real_type = typename Dune::FieldTraits<Self>::real_type;
+    using size_type = std::size_t;
+
+    enum {
+      rows = std::tuple_size<Super>::value,
+      cols = Math::max(Rows::dimension...)
+    };
+
+    template <class... Rows_,
+      REQUIRES( Concepts::Similar<Types<Rows...>, Types<Rows_...>> )>
+    MultiTypeMatrix(Rows_&&... rows)
+      : Super(std::forward<Rows_>(rows)...)
+    {}
+
+    /// Default construction of tuple of FieldVectors
+    MultiTypeMatrix() = default;
+
+    /// Construct tuple by initializing all tuple elements with a constant value
+    MultiTypeMatrix(real_type value)
+    {
+      *this = value;
+    }
+
+    /// Assignment of real number to all tuple elements
+    MultiTypeMatrix& operator=(real_type value)
+    {
+      forEach(*this, [value](auto& fv) { fv = value; });
+      return *this;
+    }
+
+    // Compound assignment operator +=
+    MultiTypeMatrix& operator+=(MultiTypeMatrix const& that)
+    {
+      forEach(range_<0,rows>, [&that,this](auto const _i) { (*this)[_i] += that[_i]; });
+      return *this;
+    }
+
+    // Compound assignment operator -=
+    MultiTypeMatrix& operator-=(MultiTypeMatrix const& that)
+    {
+      forEach(range_<0,rows>, [&that,this](auto const _i) { (*this)[_i] -= that[_i]; });
+      return *this;
+    }
+
+    // Scaling of all tuple elements by a constant value
+    MultiTypeMatrix& operator*=(real_type value)
+    {
+      forEach(*this, [value](auto& fv) { fv *= value; });
+      return *this;
+    }
+
+    // Scaling of all tuple elements by the inverse of a constant value
+    MultiTypeMatrix& operator/=(real_type value)
+    {
+      forEach(*this, [value](auto& fv) { fv /= value; });
+      return *this;
+    }
+
+    /// Const access to the tuple elements
+    template <std::size_t I, std::size_t J>
+    decltype(auto) operator()(const index_t<I>&, const index_t<J>&) const
+    {
+      return std::get<J>(std::get<I>(*this));
+    }
+
+    /// Mutable access to the tuple elements
+    template <std::size_t I, std::size_t J>
+    decltype(auto) operator()(const index_t<I>&, const index_t<J>&)
+    {
+      return std::get<J>(std::get<I>(*this));
+    }
+
+    /// Return number of elements of the tuple
+    static constexpr std::size_t num_rows()
+    {
+      return rows;
+    }
+
+    /// Return number of elements of the tuple
+    static constexpr std::size_t num_cols()
+    {
+      return cols;
+    }
+  };
+
+} // end namespace AMDiS

dune/amdis/common/MultiTypeFieldVector.hpp → dune/amdis/common/MultiTypeVector.hpp

@@ -2,93 +2,95 @@
 
 #include <tuple>
 
-#include <dune/common/ftraits.hh>
+#include <dune/functions/common/indexaccess.hh>
 
 #include <dune/amdis/common/Concepts.hpp>
+#include <dune/amdis/common/FieldTraits.hpp>
+#include <dune/amdis/common/Loops.hpp>
 #include <dune/amdis/common/Mpl.hpp>
 #include <dune/amdis/common/Size.hpp>
 
 namespace AMDiS
 {
   // forward declaration
-  template <class FV0, class... FV>
-  class MultiTypeFieldVector;
+  template <class... FV>
+  class MultiTypeVector;
 }
 
 namespace Dune
 {
-  template <class FV0, class... FV>
-  struct FieldTraits<AMDiS::MultiTypeFieldVector<FV0,FV...>>
+  template <class... FV>
+  struct FieldTraits<AMDiS::MultiTypeVector<FV...>>
   {
-    using field_type = typename Dune::FieldTraits<FV0>::field_type;
-    using real_type = typename Dune::FieldTraits<FV0>::real_type;
+    using field_type = typename AMDiS::CommonFieldTraits<FV...>::field_type;
+    using real_type = typename AMDiS::CommonFieldTraits<FV...>::real_type;
   };
 }
 
 
 namespace AMDiS
 {
-  template <class FV0, class... FV>
-  class MultiTypeFieldVector
-      : public std::tuple<FV0,FV...>
+  template <class... FV>
+  class MultiTypeVector
+      : public std::tuple<FV...>
   {
-    using Self = MultiTypeFieldVector;
-    using Super = std::tuple<FV0,FV...>;
+    using Self = MultiTypeVector;
+    using Super = std::tuple<FV...>;
 
   public:
-    using field_type = typename Dune::FieldTraits<FV0>::field_type;
-    using real_type = typename Dune::FieldTraits<FV0>::real_type;
+    using field_type = typename Dune::FieldTraits<Self>::field_type;
+    using real_type = typename Dune::FieldTraits<Self>::real_type;
     using size_type = std::size_t;
 
     enum {
       dimension = std::tuple_size<Super>::value
     };
 
-    template <class FV0_, class... FV_,
-      REQUIRES( Concepts::Similar<Types<FV0,FV...>, Types<FV0_,FV_...>> )>
-    MultiTypeFieldVector(FV0_&& fv0, FV_&&... fv)
-      : Super(std::forward<FV0_>(fv0), std::forward<FV_>(fv)...)
+    template <class... FV_,
+      REQUIRES( Concepts::Similar<Types<FV...>, Types<FV_...>> )>
+    MultiTypeVector(FV_&&... fv)
+      : Super(std::forward<FV_>(fv)...)
     {}
 
     /// Default construction of tuple of FieldVectors
-    MultiTypeFieldVector() = default;
+    MultiTypeVector() = default;
 
     /// Construct tuple by initializing all tuple elements with a constant value
-    MultiTypeFieldVector(real_type value)
+    explicit MultiTypeVector(real_type value)
     {
       *this = value;
     }
 
     /// Assignment of real number to all tuple elements
-    MultiTypeFieldVector& operator=(real_type value)
+    MultiTypeVector& operator=(real_type value)
     {
       forEach(*this, [value](auto& fv) { fv = value; });
       return *this;
     }
 
     // Compound assignment operator +=
-    MultiTypeFieldVector& operator+=(MultiTypeFieldVector const& that)
+    MultiTypeVector& operator+=(MultiTypeVector const& that)
     {
       forEach(range_<0,dimension>, [&that,this](auto const _i) { (*this)[_i] += that[_i]; });
       return *this;
     }
 
     // Compound assignment operator -=
-    MultiTypeFieldVector& operator-=(MultiTypeFieldVector const& that)
+    MultiTypeVector& operator-=(MultiTypeVector const& that)
     {
       forEach(range_<0,dimension>, [&that,this](auto const _i) { (*this)[_i] -= that[_i]; });
       return *this;
     }
 
     // Scaling of all tuple elements by a constant value
-    MultiTypeFieldVector& operator*=(real_type value)
+    MultiTypeVector& operator*=(real_type value)
     {
       forEach(*this, [value](auto& fv) { fv *= value; });
       return *this;
     }
 
     // Scaling of all tuple elements by the inverse of a constant value
-    MultiTypeFieldVector& operator/=(real_type value)
+    MultiTypeVector& operator/=(real_type value)
     {
       forEach(*this, [value](auto& fv) { fv /= value; });
       return *this;
@@ -108,36 +110,27 @@ namespace AMDiS
       return std::get<I>(*this);
     }
 
-    /// Return number of elements of the tuple
-    static constexpr std::size_t size()
+    /// Const access to the vector using multi-indices
+    template <class Index,
+      REQUIRES( Concepts::MultiIndex<Index> )>
+    decltype(auto) operator[](Index const& index) const
     {
-      return std::tuple_size<Super>::value;
+      return Dune::Functions::hybridMultiIndexAccess<field_type const&>(*this, index);
     }
 
-    // /// Return number of elements in a flat vector
-    // static constexpr std::size_t flat_size()
-    // {
-    //   return Size<Self>;
-    // }
-
-    // /// Create flat FieldVector
-    // operator Dune::FieldVector<field_type, flat_size()>() const
-    // {
-    //   Dune::FieldVector<field_type, flat_size()> result;
-
-    //   int shift = 0;
-    //   forEach(*this, [&result,&shift](auto& fv)
-    //   {
-    //     static const int sub_size = Size<decltype(fv)>;
-    //     Dune::FieldVector<field_type, sub_size> sub = fv;
-    //     for (int i = 0; i < sub_size; ++i)
-    //       result[shift + i] = sub[i];
-    //     shift += sub_size;
-    //   });
-
-    //   return result;
-    // }
+    /// Mutable access to the vector using multi-indices
+    template <class Index,
+      REQUIRES( Concepts::MultiIndex<Index> )>
+    decltype(auto) operator[](Index const& index)
+    {
+      return Dune::Functions::hybridMultiIndexAccess<field_type&>(*this, index);
+    }
 
+    /// Return number of elements of the tuple
+    static constexpr std::size_t size()
+    {
+      return dimension;
+    }
   };
 
-} // end namesüace AMDiS
+} // end namespace AMDiS

dune/amdis/common/TypeDefs.hpp

@@ -2,9 +2,7 @@
 
 #include <type_traits>
 
-#include <boost/numeric/mtl/matrix/compressed2D.hpp>
-#include <boost/numeric/mtl/vector/dense_vector.hpp>
-#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <dune/amdis/linear_algebra/Mtl.hpp>
 
 namespace AMDiS
 {

dune/amdis/common/Utility.hpp

@@ -25,6 +25,25 @@ namespace AMDiS
   using Underlying_t = typename Impl::UnderlyingType<T>::type;
 
 
+  namespace Impl
+  {
+    template <class T0, class... Ts>
+    struct CommonType
+    {
+      using type = std::common_type_t<T0, typename CommonType<Ts...>::type>;
+    };
+
+    template <class T0>
+    struct CommonType<T0>
+    {
+      using type = T0;
+    };
+  }
+
+  template <class... T>
+  using CommonType_t = typename Impl::CommonType<T...>::type;
+
+
   // ---------------------------------------------------------------------------
 
 

dune/amdis/io/FileWriterInterface.hpp

@@ -22,9 +22,6 @@ namespace AMDiS
       Parameters::get(base + "->output directory", dir_);
       name_ = "solution";
       Parameters::get(base + "->name", name_);
-
-      if (!filesystem::exists(dir_))
-        error_exit("Output directory '",dir_,"' does not exist!");
     }
 
 

dune/amdis/linear_algebra/DOFVectorBase.hpp

+#pragma once
+
+namespace AMDiS
+{
+  template <class GlobalBasisType>
+  class DOFVectorBase
+  {
+  public:
+    using GlobalBasis = GlobalBasisType;
+    using SizeType = std::size_t;
+
+    DOFVectorBase(GlobalBasisType const& basis)