From 97813675e10e9ff8a98cf33e7f43a2bca1adde08 Mon Sep 17 00:00:00 2001
From: "Praetorius, Simon" <simon.praetorius@tu-dresden.de>
Date: Tue, 3 Sep 2019 13:04:00 +0200
Subject: [PATCH] Added flat matrix and vector used for ElementMatrix and
 ElementVector

---
 src/amdis/common/CMakeLists.txt           |   2 +
 src/amdis/common/FlatMatrix.hpp           | 110 ++++++++++++++++++++++
 src/amdis/common/FlatVector.hpp           |  23 +++++
 src/amdis/linearalgebra/DOFMatrixBase.hpp |   3 +-
 src/amdis/linearalgebra/DOFVectorBase.hpp |   3 +-
 test/CMakeLists.txt                       |   3 +
 test/FlatMatVecTest.cpp                   |  96 +++++++++++++++++++
 7 files changed, 238 insertions(+), 2 deletions(-)
 create mode 100644 src/amdis/common/FlatMatrix.hpp
 create mode 100644 src/amdis/common/FlatVector.hpp
 create mode 100644 test/FlatMatVecTest.cpp

diff --git a/src/amdis/common/CMakeLists.txt b/src/amdis/common/CMakeLists.txt
index 8f5c7dd3..006a5bc0 100644
--- a/src/amdis/common/CMakeLists.txt
+++ b/src/amdis/common/CMakeLists.txt
@@ -15,6 +15,8 @@ install(FILES
     FieldMatVec.hpp
     FieldMatVec.inc.hpp
     Filesystem.hpp
+    FlatMatrix.hpp
+    FlatVector.hpp
     ForEach.hpp
     HybridSize.hpp
     Index.hpp
diff --git a/src/amdis/common/FlatMatrix.hpp b/src/amdis/common/FlatMatrix.hpp
new file mode 100644
index 00000000..7d00dc88
--- /dev/null
+++ b/src/amdis/common/FlatMatrix.hpp
@@ -0,0 +1,110 @@
+#pragma once
+
+#include <type_traits>
+#include <vector>
+
+namespace AMDiS
+{
+  /// Dense matrix with row-wise storage in flat data vector
+  template <class T, class Allocator = std::allocator<T>>
+  class FlatMatrix
+      : public std::vector<T,Allocator>
+  {
+    using Super = std::vector<T,Allocator>;
+
+  public:
+    using value_type = T;
+    using size_type = typename Super::size_type;
+    using reference = typename Super::reference;
+    using const_reference = typename Super::const_reference;
+
+  private:
+    /// Accessor to a row of the matrix.
+    template <class Vector>
+    struct AccessProxy
+    {
+      /// Access the column entry `col` of the row `row_` of this accessor.
+      /// @{
+      template <class V = Vector, std::enable_if_t<not std::is_const<V>::value,int> = 0>
+      reference operator[](size_type col)
+      {
+        return (*vector_)[row_*cols_ + col];
+      }
+
+      const_reference operator[](size_type col) const
+      {
+        return (*vector_)[row_*cols_ + col];
+      }
+      // @}
+
+      Vector* vector_;
+      size_type row_;
+      size_type cols_;
+    };
+
+  public:
+    /// Default constructor, creates an empty 0x0 matrix.
+    FlatMatrix() = default;
+
+    /// Construct a new matrix with rows `r` and columns `c` and all values
+    /// initialized by copies of value `v`.
+    FlatMatrix(size_type r, size_type c, value_type v = {})
+      : Super(r*c, v)
+      , rows_(r)
+      , cols_(c)
+    {}
+
+    /// Assign value `s` to all entries of the matrix
+    FlatMatrix& operator=(value_type s)
+    {
+      Super::assign(Super::size(), s);
+      return *this;
+    }
+
+    /// Resizes the container to contain r x c elements.
+    /**
+     * If the current `rows*cols` is greater than `r*c`, the container is reduced
+     * to its first `r*c` elements. If the current `rows*cols` is less than `r*c`,
+     * additional copies of value `v` are appended.
+     **/
+    void resize(size_type r, size_type c, value_type v = {})
+    {
+      Super::resize(r*c, v);
+      rows_ = r;
+      cols_ = c;
+    }
+
+    /// Return the number of rows of the matrix
+    size_type rows() const noexcept
+    {
+      return rows_;
+    }
+
+    /// Return the number of columns of the matrix
+    size_type cols() const noexcept
+    {
+      return cols_;
+    }
+
+    /// Return the number of entries in the matrix, i.e. `rows*cols`.
+    using Super::size;
+
+    /// Return accessor to the `row` of the matrix, see \ref AccessProxy
+    /// @{
+    AccessProxy<Super> operator[](size_type row)
+    {
+      return AccessProxy<Super>{static_cast<Super*>(this), row, cols_};
+    }
+
+    AccessProxy<Super const> operator[](size_type row) const
+    {
+      return AccessProxy<Super const>{static_cast<Super const*>(this), row, cols_};
+    }
+    /// @}
+
+  private:
+    size_type rows_ = 0;
+    size_type cols_ = 0;
+  };
+
+} // end namespace AMDiS
diff --git a/src/amdis/common/FlatVector.hpp b/src/amdis/common/FlatVector.hpp
new file mode 100644
index 00000000..079cbae6
--- /dev/null
+++ b/src/amdis/common/FlatVector.hpp
@@ -0,0 +1,23 @@
+#pragma once
+
+#include <vector>
+
+namespace AMDiS
+{
+  /// Flat data vector to be used in assembling as element vector
+  template <class T, class Allocator = std::allocator<T>>
+  class FlatVector
+      : public std::vector<T,Allocator>
+  {
+  public:
+    using std::vector<T,Allocator>::vector;
+
+    /// Assign value `s` to all entries of the vector
+    FlatVector& operator=(T s)
+    {
+      this->assign(this->size(), s);
+      return *this;
+    }
+  };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/DOFMatrixBase.hpp b/src/amdis/linearalgebra/DOFMatrixBase.hpp
index a197930c..181069c6 100644
--- a/src/amdis/linearalgebra/DOFMatrixBase.hpp
+++ b/src/amdis/linearalgebra/DOFMatrixBase.hpp
@@ -5,6 +5,7 @@
 #include <dune/common/timer.hh>
 
 #include <amdis/OperatorList.hpp>
+#include <amdis/common/FlatMatrix.hpp>
 #include <amdis/common/Math.hpp>
 #include <amdis/typetree/MultiIndex.hpp>
 #include <amdis/typetree/TreePath.hpp>
@@ -43,7 +44,7 @@ namespace AMDiS
     using BaseMatrix = typename Backend::BaseMatrix;
 
     /// The type of the matrix filled on an element with local contributions
-    using ElementMatrix = Dune::DynamicMatrix<value_type>;
+    using ElementMatrix = FlatMatrix<value_type>;
 
   public:
     /// Constructor. Stores the shared_ptr to the bases.
diff --git a/src/amdis/linearalgebra/DOFVectorBase.hpp b/src/amdis/linearalgebra/DOFVectorBase.hpp
index 50f6d6a7..ad33a91a 100644
--- a/src/amdis/linearalgebra/DOFVectorBase.hpp
+++ b/src/amdis/linearalgebra/DOFVectorBase.hpp
@@ -9,6 +9,7 @@
 #include <amdis/DataTransfer.hpp>
 #include <amdis/GridTransferManager.hpp>
 #include <amdis/OperatorList.hpp>
+#include <amdis/common/FlatVector.hpp>
 #include <amdis/common/Math.hpp>
 #include <amdis/common/TypeTraits.hpp>
 #include <amdis/linearalgebra/DOFVectorInterface.hpp>
@@ -73,7 +74,7 @@ namespace AMDiS
     using BaseVector = typename Backend::BaseVector;
 
     /// The type of the vector filled on an element with local contributions
-    using ElementVector = Dune::DynamicVector<value_type>;
+    using ElementVector = FlatVector<value_type>;
 
     /// Defines an interface to transfer the data during grid adaption
     using DataTransfer = DataTransferInterface<Self>;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 3d6c4220..d5e76ab8 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -45,6 +45,9 @@ dune_add_test(SOURCES FiniteElementTypeTest.cpp
 dune_add_test(SOURCES FilesystemTest.cpp
   LINK_LIBRARIES amdis)
 
+dune_add_test(SOURCES FlatMatVecTest.cpp
+  LINK_LIBRARIES amdis)
+
 dune_add_test(SOURCES GlobalIdSetTest.cpp
   LINK_LIBRARIES amdis
   MPI_RANKS 2
diff --git a/test/FlatMatVecTest.cpp b/test/FlatMatVecTest.cpp
new file mode 100644
index 00000000..b02816a3
--- /dev/null
+++ b/test/FlatMatVecTest.cpp
@@ -0,0 +1,96 @@
+#include <cmath>
+
+#include <amdis/common/FlatMatrix.hpp>
+#include <amdis/common/FlatVector.hpp>
+
+#include "Tests.hpp"
+
+using namespace AMDiS;
+
+// -----------------------------------------------------------------------------
+void test0()
+{
+  using V = FlatVector<double>;
+  // create an empty vector
+  V a;
+
+  // size of the vector
+  AMDIS_TEST_EQ(a.size(), 0);
+
+  // Resize the vector
+  a.resize(3, 0.0);
+  AMDIS_TEST_EQ(a.size(), 3);
+  AMDIS_TEST_EQ(a[0], 0.0);
+  AMDIS_TEST_EQ(a[1], 0.0);
+  AMDIS_TEST_EQ(a[2], 0.0);
+
+  // access component
+  a[0] = 0.0;
+  a[1] = 0.5;
+  a[2] = 1.0;
+
+  // assign value to all entries
+  a = 2.0;
+  AMDIS_TEST_EQ(a[0], 2.0);
+  AMDIS_TEST_EQ(a[1], 2.0);
+  AMDIS_TEST_EQ(a[2], 2.0);
+}
+
+// -----------------------------------------------------------------------------
+void test1()
+{
+  using M = FlatMatrix<double>;
+  // create an empty matrix
+  M m;
+  AMDIS_TEST_EQ(m.size(), 0);
+  AMDIS_TEST_EQ(m.rows(), 0);
+  AMDIS_TEST_EQ(m.cols(), 0);
+
+  // resize the matrix
+  m.resize(2,3, 0.0);
+  AMDIS_TEST_EQ(m.size(), 2*3);
+  AMDIS_TEST_EQ(m.rows(), 2);
+  AMDIS_TEST_EQ(m.cols(), 3);
+
+  typename M::difference_type s = std::distance(m.begin(),m.end());
+  AMDIS_TEST_EQ(s, 2*3);
+
+  // test element access
+  AMDIS_TEST_EQ(m[0][0], 0.0);
+  AMDIS_TEST_EQ(m[1][2], 0.0);
+
+  // test assignment of value
+  m = 1.0;
+  AMDIS_TEST_EQ(m[0][0], 1.0);
+  AMDIS_TEST_EQ(m[1][2], 1.0);
+
+  for (std::size_t i = 0; i < m.rows(); ++i)
+    for (std::size_t j = 0; j < m.cols(); ++j)
+      m[i][j] = 2.0;
+
+  // access matrix entries
+  for (std::size_t i = 0; i < m.rows(); ++i)
+    for (std::size_t j = 0; j < m.cols(); ++j)
+      AMDIS_TEST_EQ(m[i][j], 2.0);
+
+  // test iterators
+  for (auto& mij : m)
+    mij = 3.0;
+
+  for (auto const& mij : m)
+    AMDIS_TEST_EQ(mij, 3.0);
+
+  // test data array
+  m = 4.0;
+  double const* data = m.data();
+  for (std::size_t i = 0; i < m.size(); ++i)
+    AMDIS_TEST_EQ(data[i], 4.0);
+}
+
+int main()
+{
+  test0();
+  test1();
+
+  return report_errors();
+}
-- 
GitLab