Merge pull request #1312 from su2code/gradient_type

Another charge against pointer to pointer

Author: pcarruscag

Common/include/containers/container_decorators.hpp

@@ -30,6 +30,43 @@
 
 #include "C2DContainer.hpp"
 
+/*!
+ * \brief Class to represent a matrix (without owning the data, this just wraps a pointer).
+ */
+template<class T>
+class CMatrixView {
+public:
+  using Scalar = typename std::remove_const<T>::type;
+  using Index = unsigned long;
+
+private:
+  T* m_ptr;
+  Index m_cols;
+
+public:
+  CMatrixView(T* ptr = nullptr, Index cols = 0) : m_ptr(ptr), m_cols(cols) {}
+
+  template<class U> friend class CMatrixView;
+  template<class U>
+  CMatrixView(const CMatrixView<U>& other) : m_ptr(other.m_ptr), m_cols(other.m_cols) {}
+
+  explicit CMatrixView(su2matrix<Scalar>& mat) : m_ptr(mat.data()), m_cols(mat.cols()) {}
+
+  template<class U = T, su2enable_if<std::is_const<U>::value> = 0>
+  explicit CMatrixView(const su2matrix<Scalar>& mat) : m_ptr(mat.data()), m_cols(mat.cols()) {}
+
+  const Scalar* operator[] (Index i) const noexcept { return &m_ptr[i*m_cols]; }
+  const Scalar& operator() (Index i, Index j) const noexcept { return m_ptr[i*m_cols + j]; }
+
+  template<class U = T, su2enable_if<!std::is_const<U>::value> = 0>
+  Scalar* operator[] (Index i) noexcept { return &m_ptr[i*m_cols]; }
+
+  template<class U = T, su2enable_if<!std::is_const<U>::value> = 0>
+  Scalar& operator() (Index i, Index j) noexcept { return m_ptr[i*m_cols + j]; }
+
+  friend CMatrixView operator+ (CMatrixView mv, Index incr) { return CMatrixView(mv[incr], mv.m_cols); }
+};
+
 /*!
  * \class C3DContainerDecorator
  * \brief Decorate a matrix type (Storage) with 3 dimensions.
@@ -44,6 +81,9 @@ class C3DContainerDecorator {
   static constexpr bool IsRowMajor = true;
   static constexpr bool IsColumnMajor = false;
 
+  using Matrix = CMatrixView<Scalar>;
+  using ConstMatrix = CMatrixView<const Scalar>;
+
   using CInnerIter = typename Storage::CInnerIter;
   template<class T, size_t N>
   using CInnerIterGather = typename Storage::template CInnerIterGather<simd::Array<T,N> >;
@@ -78,6 +118,18 @@ class C3DContainerDecorator {
   Scalar& operator() (Index i, Index j, Index k) noexcept { return m_storage(i, j*m_innerSz + k); }
   const Scalar& operator() (Index i, Index j, Index k) const noexcept { return m_storage(i, j*m_innerSz + k); }
 
+  /*!
+   * \brief Matrix access.
+   */
+  Matrix operator[] (Index i) noexcept { return Matrix(m_storage[i], m_innerSz); }
+  ConstMatrix operator[] (Index i) const noexcept { return ConstMatrix(m_storage[i], m_innerSz); }
+
+  /*!
+   * \brief Matrix access with an offset.
+   */
+  Matrix operator() (Index i, Index j) noexcept { return Matrix(m_storage[i]+j*m_innerSz, m_innerSz); }
+  ConstMatrix operator() (Index i, Index j) const noexcept { return ConstMatrix(m_storage[i]+j*m_innerSz, m_innerSz); }
+
   /*!
    * \brief Get a scalar iterator to the inner-most dimension of the container.
    */
@@ -105,42 +157,11 @@ class C3DContainerDecorator {
 };
 
 /*!
- * \brief Some typedefs for the
+ * \brief Some typedefs for 3D containers
  */
 using C3DIntMatrix = C3DContainerDecorator<su2matrix<unsigned long> >;
 using C3DDoubleMatrix = C3DContainerDecorator<su2activematrix>;
-
-/*!
- * \class CVectorOfMatrix
- * \brief This contrived container is used to store small matrices in a contiguous manner
- *        but still present the "su2double**" interface to the outside world.
- *        The "interface" part should be replaced by something more efficient, e.g. a "matrix view".
- */
-class CVectorOfMatrix: public C3DDoubleMatrix {
-private:
-  su2matrix<Scalar*> interface;
-
-public:
-  CVectorOfMatrix() = default;
-
-  CVectorOfMatrix(Index length, Index rows, Index cols, Scalar value = 0) noexcept {
-    resize(length, rows, cols, value);
-  }
-
-  void resize(Index length, Index rows, Index cols, Scalar value = 0) noexcept {
-    C3DDoubleMatrix::resize(length, rows, cols, value);
-    interface.resize(length,rows);
-    for(Index i=0; i<length; ++i)
-      for(Index j=0; j<rows; ++j)
-        interface(i,j) = &(*this)(i,j,0);
-  }
-
-  /*!
-   * \brief Matrix-wise access.
-   */
-  Scalar** operator[] (Index i) noexcept { return interface[i]; }
-  const Scalar* const* operator[] (Index i) const noexcept { return interface[i]; }
-};
+using CVectorOfMatrix = C3DDoubleMatrix;
 
 /*!
  * \class C2DDummyLastView

Common/include/geometry/dual_grid/CPoint.hpp

@@ -809,7 +809,7 @@ class CPoint {
    * \param[in] iPoint - Index of the point.
    * \return Grid velocity gradient at the point.
    */
-  inline const su2double* const* GetGridVel_Grad(unsigned long iPoint) const { return GridVel_Grad[iPoint]; }
+  inline CMatrixView<const su2double> GetGridVel_Grad(unsigned long iPoint) const { return GridVel_Grad[iPoint]; }
 
   /*!
    * \brief Set the adjoint values of the (geometric) coordinates.