added more tests and corrected MultiTypeMatrix

dune/amdis/common/Mpl.hpp

@@ -153,8 +153,21 @@ namespace AMDiS
   constexpr bool_t<!B> operator!(bool_t<B>) { return {}; }
 
 
-  template <class T, T value0, T... values>
-  using is_equal = all_of_t<(value0 == values)...>;
+  namespace Impl
+  {
+    template <class T, T I0, T... Is>
+    struct IsEqual
+    {
+      template <T I>
+      using compare = bool_t<I0 == I>;
+
+      enum { value = all_of_v<compare<Is>::value...> };
+    };
+
+  } // end namespace Impl
+
+  template <class T, T... values>
+  using is_equal = Impl::IsEqual<T,values...>;
 
   template <class T, class... Ts>
   using is_one_of = or_t<std::is_same<T, Ts>::value...>;

dune/amdis/common/MultiTypeMatrix.hpp

@@ -46,10 +46,8 @@ namespace AMDiS
     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...)
-    };
+    static constexpr int rows = std::tuple_size<Super>::value;
+    static constexpr int cols = Math::max(Rows::dimension...);
 
     template <class... Rows_,
       REQUIRES( Concepts::Similar<Types<Rows...>, Types<Rows_...>> )>
@@ -115,6 +113,21 @@ namespace AMDiS
       return std::get<J>(std::get<I>(*this));
     }
 
+
+    /// Const access to the rows
+    template <std::size_t I>
+    decltype(auto) operator[](const index_t<I>&) const
+    {
+      return std::get<I>(*this);
+    }
+
+    /// Mutable access to the rows
+    template <std::size_t I>
+    decltype(auto) operator[](const index_t<I>&)
+    {
+      return std::get<I>(*this);
+    }
+
     /// Return number of elements of the tuple
     static constexpr std::size_t num_rows()
     {
@@ -126,6 +139,11 @@ namespace AMDiS
     {
       return cols;
     }
+
+    static constexpr std::size_t size()
+    {
+      return rows;
+    }
   };
 
 } // end namespace AMDiS

dune/amdis/common/MultiTypeVector.hpp

@@ -42,9 +42,7 @@ namespace AMDiS
     using real_type = typename Dune::FieldTraits<Self>::real_type;
     using size_type = std::size_t;
 
-    enum {
-      dimension = std::tuple_size<Super>::value
-    };
+    static constexpr int dimension = std::tuple_size<Super>::value;
 
     template <class... FV_,
       REQUIRES( Concepts::Similar<Types<FV...>, Types<FV_...>> )>

dune/amdis/common/Size.hpp

@@ -54,5 +54,4 @@ namespace AMDiS
   template <class T>
   constexpr std::size_t Size = Impl::SizeImpl<std::decay_t<T>>::value;
 
-
 } // end namespace AMDiS

dune/amdis/common/Tags.hpp

+#pragma once
+
+namespace AMDiS
+{
+  namespace tag
+  {
+    // some tags representing traits classes or categories
+    struct none {};
+    struct scalar {};
+    struct vector {};
+    struct matrix {};
+    struct unknown {};
+
+  } // end namespace tag
+
+} // end namespace AMDiS

dune/amdis/common/ValueCategory.hpp

@@ -6,20 +6,10 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 
+#include <dune/amdis/common/Tags.hpp>
+
 namespace AMDiS
 {
-  namespace tag
-  {
-    // some tags representing traits classes or categories
-    struct none {};
-    struct scalar {};
-    struct vector {};
-    struct matrix {};
-    struct unknown {};
-
-  } // end namespace tag
-
-
   /// Category of type T, e.g. scalar, vector matrix, specified by a tag
   template <class T, class = void>
   struct ValueCategory

dune/amdis/utility/FiniteElementType.hpp

 #pragma once
 
-#include <dune/amdis/common/ValueCategory.hpp>
+#include <cmath>
+
+#include <dune/amdis/common/Loops.hpp>
+#include <dune/amdis/common/Mpl.hpp>
+#include <dune/amdis/common/Tags.hpp>
 
 namespace AMDiS
 {

dune/amdis/utility/MultiIndex.hpp

@@ -14,45 +14,4 @@ namespace AMDiS
 
   std::size_t flatMultiIndex(Dune::ReservedVector<std::size_t,1> const& idx) { return idx[0]; }
 
-
-
-  template <int rank, class Index>
-  struct MultiIndexResolver
-  {
-    MultiIndexResolver(Index const& index)
-      : index_(index)
-    {}
-
-    template <class C>
-    decltype(auto) operator()(C&& c)
-    {
-      return resolve(std::forward<C>(c), int_<rank>);
-    }
-
-  private:
-    template <class C, int r>
-    decltype(auto) resolve(C&& c, int_t<r>)
-    {
-      using namespace Dune::Indices;
-      auto&& subIndex = Dune::Functions::shiftedMultiIndex(index_);
-      auto&& subIndexResolver = MultiIndexResolver<r-1, decltype(subIndex)>(subIndex);
-      return Dune::Functions::hybridIndexAccess(c, index_[_0], subIndexResolver);
-    }
-
-    template <class C>
-    decltype(auto) resolve(C&& c, int_t<0>)
-    {
-      return std::forward<C>(c);
-    }
-
-    Index const& index_;
-  };
-
-  template <int rank, class Container, class MultiIndex>
-  decltype(auto) multiIndexAccess(Container&& c, MultiIndex const& index)
-  {
-    MultiIndexResolver<rank, MultiIndex> multiIndexResolver(index);
-    return multiIndexResolver(c);
-  }
-
 } // end namespace AMDiS

dune/amdis/utility/RangeType.hpp

@@ -2,6 +2,8 @@
 
 #include <type_traits>
 
+#include <dune/common/fvector.hh>
+
 #include <dune/amdis/common/IndexSeq.hpp>
 #include <dune/amdis/common/MultiTypeVector.hpp>
 
@@ -29,8 +31,6 @@ namespace AMDiS
     {
       using LocalBasis = typename Node::FiniteElement::Traits::LocalBasisType;
       using type = typename LocalBasis::Traits::RangeType;
-      // using type_ = typename LocalBasis::Traits::RangeType;
-      // using type = std::conditional_t<std::is_same<type_, Dune::FieldVector<double,1>>::value, double, type_>;
     };
 
     // Power node

dune/amdis/utility/Visitor.hpp

 #pragma once
 
 #include <dune/typetree/visitor.hh>
+#include <dune/amdis/common/Mpl.hpp>
 #include <dune/amdis/utility/Traversal.hpp>
 
 namespace AMDiS

test/CMakeLists.txt

@@ -9,11 +9,23 @@ dune_add_test(SOURCES ExpressionsTest.cpp
 dune_add_test(SOURCES FieldMatVecTest.cpp
   LINK_LIBRARIES amdis)
 
+dune_add_test(SOURCES FiniteElementTypeTest.cpp
+  LINK_LIBRARIES amdis)
+
 dune_add_test(SOURCES FilesystemTest.cpp
   LINK_LIBRARIES amdis)
 
 dune_add_test(SOURCES MultiTypeVectorTest.cpp
   LINK_LIBRARIES amdis)
 
+dune_add_test(SOURCES MultiTypeMatrixTest.cpp
+  LINK_LIBRARIES amdis)
+
+dune_add_test(SOURCES RangeTypeTest.cpp
+  LINK_LIBRARIES amdis)
+
 dune_add_test(SOURCES StringTest.cpp
   LINK_LIBRARIES amdis)
+
+dune_add_test(SOURCES TreeDataTest.cpp
+  LINK_LIBRARIES amdis)

test/FiniteElementTypeTest.cpp

+#include <dune/common/filledarray.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
+#include <dune/amdis/utility/FiniteElementType.hpp>
+#include "Tests.hpp"
+
+int main()
+{
+  using namespace AMDiS;
+  using namespace Dune;
+  using namespace Dune::Functions;
+  using namespace Dune::Functions::BasisBuilder;
+  using namespace Dune::Indices;
+
+  // create grid
+  FieldVector<double, 2> L; L = 1.0;
+  auto s = Dune::filledArray<2>(1);
+  YaspGrid<2> grid(L, s);
+  auto gridView = grid.leafGridView();
+
+  // create basis
+  static const int k = 1;
+  auto taylorHoodBasis = makeBasis(
+    gridView,
+    composite(
+      power<2>(lagrange<k+1>()),
+      lagrange<k>()
+    ));
+
+  auto localView = taylorHoodBasis.localView();
+
+  for (auto const& e : elements(gridView)) {
+    localView.bind(e);
+
+    auto node = localView.tree();
+    AMDIS_TEST_EQ( polynomialDegree(node), k+1 ); // maximum over all polynomial degrees
+    static_assert( std::is_same<tag::unknown, FiniteElementType_t<decltype(node)>>::value, "" );
+
+    auto v_node = TypeTree::child(node, _0);
+    AMDIS_TEST_EQ( polynomialDegree(v_node), k+1 );
+
+    auto p_node = TypeTree::child(node, _1);
+    AMDIS_TEST_EQ( polynomialDegree(p_node), k );
+  }
+
+  return report_errors();
+}

test/MultiTypeMatrixTest.cpp

+#include <dune/amdis/common/FieldMatVec.hpp>
+#include <dune/amdis/common/MultiTypeMatrix.hpp>
+
+using namespace AMDiS;
+
+int main()
+{
+  using V0 = MultiTypeVector<FieldMatrix<double,2,2>, FieldMatrix<double,2,3>>;
+  using V1 = MultiTypeVector<FieldMatrix<double,3,2>, FieldMatrix<double,3,3>>;
+  using M = MultiTypeMatrix<V0, V1>;
+
+  // constructors
+  M m0{}; // default constructor
+
+  V0 v0{};
+  V1 v1{};
+
+  M m1{v0,v1}; // constructor with Vi arguments
+  M m2{m1}; // copy constructor
+  m0 = m1; // copy assignment
+
+  M m3{0.0}; // construction from constant
+  m3 = 1.0; // assignment of constant
+
+  {
+    M m0_tmp;
+    M m1_tmp{std::move(m0_tmp)}; // move constructor
+    m1 = std::move(m1_tmp); // move assignment
+  }
+
+  using namespace Dune::Indices;
+
+  // element access
+  DUNE_UNUSED auto fm0 = m1(_0,_0); // const block access
+  DUNE_UNUSED auto fm1 = m1[_1][_1]; // const block access
+  m2(_0,_0) = fm0; // mutable block access
+  m2[_1][_1] = fm1; // mutable block access
+
+  DUNE_UNUSED double x = m2(_0,_0)[1][1];  // const access
+
+  // compound assignment operators
+  m1 += m2;
+  m1 *= 2.0;
+  m1 /= 2.0;
+
+}
\ No newline at end of file

test/RangeTypeTest.cpp

+#include <dune/common/filledarray.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
+#include <dune/amdis/utility/RangeType.hpp>
+#include "Tests.hpp"
+
+int main()
+{
+  using namespace AMDiS;
+  using namespace Dune;
+  using namespace Dune::Functions;
+  using namespace Dune::Functions::BasisBuilder;
+  using namespace Dune::Indices;
+
+  // create grid
+  FieldVector<double, 2> L; L = 1.0;
+  auto s = Dune::filledArray<2>(1);
+  YaspGrid<2> grid(L, s);
+  auto gridView = grid.leafGridView();
+
+  // create basis
+  static const int k = 1;
+  auto taylorHoodBasis = makeBasis(
+    gridView,
+    composite(
+      power<2>(lagrange<k+1>()),
+      lagrange<k>()
+    ));
+
+  auto localView = taylorHoodBasis.localView();
+
+  for (auto const& e : elements(gridView)) {
+    localView.bind(e);
+
+    auto node = localView.tree();
+    using NodeRange = MultiTypeVector<FieldVector<double,2>, FieldVector<double,1>>;
+    static_assert( std::is_same<NodeRange, RangeType_t<decltype(node)>>::value, "" );
+
+    auto v_node = TypeTree::child(node, _0);
+    using VNodeRange = FieldVector<double,2>;
+    static_assert( std::is_same<VNodeRange, RangeType_t<decltype(v_node)>>::value, "" );
+
+    auto p_node = TypeTree::child(node, _1);
+    using PNodeRange = FieldVector<double,1>;
+    static_assert( std::is_same<PNodeRange, RangeType_t<decltype(p_node)>>::value, "" );
+  }
+
+  return report_errors();
+}

test/TreeDataTest.cpp

+#include <iostream>
+#include <memory>
+#include <functional>
+
+#include <dune/common/filledarray.hh>
+#include <dune/grid/yaspgrid.hh>
+
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <dune/amdis/utility/TreeData.hpp>
+
+#include "Tests.hpp"
+
+using namespace AMDiS;
+
+
+template <class Node>
+using NodeData = double;
+
+template <class Tree>
+bool operator==(TreeData<Tree,NodeData,false> const& t1, TreeData<Tree,NodeData,false> const& t2)
+{
+  AMDIS_TEST(t1.tree() == t2.tree() && t1.tree() != nullptr);
+
+  bool same = true;
+  forEachNode(*t1.tree(), [&](auto const& node, auto) {
+    same = same && (t1[node] == t2[node]);
+  });
+
+  return same;
+}
+
+template <class Tree>
+bool operator==(TreeData<Tree,NodeData,true> const& t1, TreeData<Tree,NodeData,true> const& t2)
+{
+  AMDIS_TEST(t1.tree() == t2.tree() && t1.tree() != nullptr);
+
+  bool same = true;
+  forEachLeafNode(*t1.tree(), [&](auto const& node, auto) {
+    same = same && (t1[node] == t2[node]);
+  });
+
+  return same;
+}
+
+int main ()
+{
+  using namespace Dune;
+  using namespace Dune::Functions;
+  using namespace Dune::Functions::BasisBuilder;
+
+  FieldVector<double, 2> L; L = 1.0;
+  auto s = Dune::filledArray<2>(1);
+  YaspGrid<2> grid(L, s);
+  auto gridView = grid.leafGridView();
+
+  auto taylorHoodBasis = makeBasis(
+    gridView,
+    composite(
+      power<2>(lagrange<2>()),
+      lagrange<1>()
+    ));
+
+  auto localView = taylorHoodBasis.localView();
+  auto const& tree = localView.tree();
+  using Tree = std::remove_reference_t<decltype(tree)>;
+
+  // test treeData for all nodes
+  {
+    // call default-constructor
+    TreeData<Tree, NodeData, false> treeData;
+    treeData.init(tree);
+
+    forEachNode(tree, [&](auto const& node, auto) {
+      treeData[node] = double(node.treeIndex());
+    });
+
+    // call constructor with tree
+    TreeData<Tree, NodeData, false> treeData1(tree);
+
+    // call init on non-empty treeData
+    treeData1.init(tree);
+
+    // call copy-constructor
+    TreeData<Tree, NodeData, false> treeData2(treeData);
+    AMDIS_TEST(treeData == treeData2);
+
+    // call copy-assignment operator on empty treeData
+    TreeData<Tree, NodeData, false> treeData3;
+    treeData3 = treeData;
+    AMDIS_TEST(treeData == treeData3);
+
+    // call copy-assignment operator on non-empty treeData
+    treeData2 = treeData3;
+    AMDIS_TEST(treeData3 == treeData2);
+
+    // call move-assignment operator on non-empty treeData
+    treeData = std::move(treeData2);
+
+    // call move-constructor
+    TreeData<Tree, NodeData, false> treeData4(std::move(treeData3));
+  }
+
+  // test treeData for leaf only
+  {
+    // call default-constructor
+    TreeData<Tree, NodeData, true> treeData;
+    treeData.init(tree);
+
+    forEachLeafNode(tree, [&](auto const& node, auto) {
+      treeData[node] = double(node.treeIndex());
+    });
+
+    // call constructor with tree
+    TreeData<Tree, NodeData, true> treeData1(tree);
+
+    // call init on non-empty treeData
+    treeData1.init(tree);
+
+    // call copy-constructor
+    TreeData<Tree, NodeData, true> treeData2(treeData);
+    AMDIS_TEST(treeData == treeData2);
+
+    // call copy-assignment operator on empty treeData
+    TreeData<Tree, NodeData, true> treeData3;
+    treeData3 = treeData;
+    AMDIS_TEST(treeData == treeData3);
+
+    // call copy-assignment operator on non-empty treeData
+    treeData2 = treeData3;
+    AMDIS_TEST(treeData3 == treeData2);
+
+    // call move-assignment operator on non-empty treeData
+    treeData = std::move(treeData2);
+
+    // call move-constructor
+    TreeData<Tree, NodeData, true> treeData4(std::move(treeData3));
+  }
+
+  // test for operations with uninitialized tree
+  {
+    // call default-constructor without initialization
+    TreeData<Tree, NodeData, true> treeData;
+
+    // call copy-constructor
+    TreeData<Tree, NodeData, true> treeData2(treeData);
+
+    // call move-constructor
+    TreeData<Tree, NodeData, true> treeData3(std::move(treeData));
+  }
+
+  return report_errors();
+}