Tensor3D should "have-a" vector, but not "be-a" vector.

- Tensor3D wraps a vector<T> member rather than directly inheriting from the class
- Fixes a typo in the field tests when converting old -> new syntax for accessing Tensor3Ds

tdms/include/arrays/tensor3d.h

@@ -22,18 +22,31 @@
  * @tparam T Numerical datatype
  */
 template<typename T>
-class Tensor3D : public std::vector<T> {
+class Tensor3D {
+private:
+  /** @brief Convert a 3D (i,j,k) index to the corresponding index in the
+   * strided storage. */
+  int to_global_index(int i, int j, int k) const {
+    return i * n_layers_ * n_cols_ + j * n_layers_ + k;
+  }
+  int to_global_index(const ijk &index_3d) const {
+    return to_global_index(index_3d.i, index_3d.j, index_3d.k);
+  }
+
 protected:
   int n_layers_ = 0;
   int n_cols_ = 0;
   int n_rows_ = 0;
 
+  /*! Strided vector that will store the array data */
+  std::vector<T> data_;
+
   /*! The total number of elements, according to the dimension values currently
    * set. */
   int total_elements() const { return n_layers_ * n_cols_ * n_rows_; }
 
 public:
-  Tensor3D() : std::vector<T>(){};
+  Tensor3D() = default;
   Tensor3D(int n_layers, int n_cols, int n_rows) {
     allocate(n_layers, n_cols, n_rows);
   }
@@ -43,28 +56,26 @@ class Tensor3D : public std::vector<T> {
 
   /** @brief Subscript operator for the Tensor, retrieving the (i,j,k)-th
    * element. */
-  T &operator()(int i, int j, int k) {
-    return *(this->begin() + i * n_layers_ * n_cols_ + j * n_layers_ + k);
-  }
+  T &operator()(int i, int j, int k) { return data_[to_global_index(i, j, k)]; }
   T operator()(int i, int j, int k) const {
-    return *(this->begin() + i * n_layers_ * n_cols_ + j * n_layers_ + k);
+    return data_[to_global_index(i, j, k)];
   }
 
   /** @brief Subscript operator for the Tensor, retrieving the (i,j,k)-th
    * element. */
   T &operator[](const ijk &index_3d) {
-    return this->operator()(index_3d.i, index_3d.j, index_3d.k);
+    return data_[to_global_index(index_3d)];
   }
   T operator[](const ijk &index_3d) const {
-    return this->operator()(index_3d.i, index_3d.j, index_3d.k);
+    return data_[to_global_index(index_3d)];
   }
 
   /**
    * @brief Whether the tensor contains any elements.
    * @details Strictly speaking, checks whether the total_elements() method
    * returns 0, indicating that this tensor has no size and thus no elements.
    */
-  bool has_elements() const { return this->total_elements() != 0; }
+  bool has_elements() const { return total_elements() != 0; }
 
   /**
    * @brief Allocate memory for this tensor given the dimensions passed.
@@ -76,7 +87,7 @@ class Tensor3D : public std::vector<T> {
     n_layers_ = n_layers;
     n_cols_ = n_cols;
     n_rows_ = n_rows;
-    this->resize(this->total_elements());
+    data_.resize(total_elements());
   }
 
   /**
@@ -93,28 +104,28 @@ class Tensor3D : public std::vector<T> {
    */
   void initialise(T ***buffer, int n_layers, int n_cols, int n_rows,
                   bool buffer_leads_n_layers = false) {
-    this->allocate(n_layers, n_cols, n_rows);
+    allocate(n_layers, n_cols, n_rows);
     if (buffer_leads_n_layers) {
       for (int k = 0; k < n_layers_; k++) {
         for (int j = 0; j < n_cols_; j++) {
           for (int i = 0; i < n_rows_; i++) {
-            this->operator()(i, j, k) = buffer[k][j][i];
+            operator()(i, j, k) = buffer[k][j][i];
           }
         }
       }
     } else {
       for (int k = 0; k < n_layers_; k++) {
         for (int j = 0; j < n_cols_; j++) {
           for (int i = 0; i < n_rows_; i++) {
-            this->operator()(i, j, k) = buffer[i][j][k];
+            operator()(i, j, k) = buffer[i][j][k];
           }
         }
       }
     }
   }
 
   /** @brief Set all elements in the tensor to 0. */
-  void zero() { std::fill(this->begin(), this->end(), 0); }
+  void zero() { std::fill(data_.begin(), data_.end(), 0); }
 
   /**
    * @brief Computes the Frobenius norm of the tensor.
@@ -130,7 +141,7 @@ class Tensor3D : public std::vector<T> {
    */
   double frobenius() const {
     T norm_val = 0;
-    for (const T &element_value : *this) {
+    for (const T &element_value : data_) {
       norm_val += std::abs(element_value) * std::abs(element_value);
     }
     return std::sqrt(norm_val);

tdms/tests/unit/field_tests/test_Field_interpolation.cpp

@@ -26,12 +26,12 @@ using tdms_tests::TOLERANCE;
 We will test the performance of BLi at interpolating a known field to the centre
 of each Yee cell, for both the E- and H- fields. The benchmark for success will
 be superior or equivalent error (Frobenius and 1D-dimension-wise norm) to that
-produced by MATLAB performing the same functionality. Frobenious error is the
+produced by MATLAB performing the same functionality. Frobenius error is the
 Frobenius norm of the 3D matrix whose (i,j,k)-th entry is the error in the
 interpolated value at Yee cell centre (i,j,k). The slice error, for a fixed j,k,
 is the norm-error of the 1D array of interpolated values along the axis (:,j,k).
 The maximum of this is then the max-slice error: keeping track of this ensures
-us that the behaviour of BLi is consistent, and does not dramaically
+us that the behaviour of BLi is consistent, and does not dramatically
 over-compensate in some areas and under-compensate in others.
 
 All tests will be performed with cell sizes Dx = 0.25, Dy = 0.1, Dz = 0.05, over
@@ -77,7 +77,7 @@ TEST_CASE("E-field interpolation check") {
 
   // setup the "split" E-field components
   ElectricSplitField E_split(Nx - 1, Ny - 1, Nz - 1);
-  E_split.allocate();// alocates Nx, Ny, Nz memory space here
+  E_split.allocate();// allocates Nx, Ny, Nz memory space here
   E_split.tot +=
           1;// correct the "number of datapoints" variable for these fields
   // setup for non-split field components
@@ -127,14 +127,12 @@ TEST_CASE("E-field interpolation check") {
   Ez_exact[k][j][i] is the field component at position (x_lower,y_lower,z_lower)
   + (i+0.5,j+0.5,k+0.5)*cellDims
   */
-  Tensor3D<double> Ex_error, Ex_split_error, Ey_error, Ey_split_error, Ez_error,
-          Ez_split_error;
-  Ex_error.allocate(Nz, Ny, Nx - 1);
-  Ex_split_error.allocate(Nz, Ny, Nx - 1);
-  Ey_error.allocate(Nz, Ny - 1, Nx);
-  Ey_split_error.allocate(Nz, Ny - 1, Nx);
-  Ez_error.allocate(Nz - 1, Ny, Nx);
-  Ez_split_error.allocate(Nz - 1, Ny, Nx);
+  Tensor3D<double> Ex_error(Nz, Ny, Nx - 1);
+  Tensor3D<double> Ex_split_error(Nz, Ny, Nx - 1);
+  Tensor3D<double> Ey_error(Nz, Ny - 1, Nx);
+  Tensor3D<double> Ey_split_error(Nz, Ny - 1, Nx);
+  Tensor3D<double> Ez_error(Nz - 1, Ny, Nx);
+  Tensor3D<double> Ez_split_error(Nz - 1, Ny, Nx);
   // now interpolate
   // note that we aren't interpolating to position 0 (before 1st point) or
   // N{x,y,z} (after last point)
@@ -180,8 +178,8 @@ TEST_CASE("E-field interpolation check") {
           double Ez_interp =
                   E.interpolate_to_centre_of(AxialDirection::Z, current_cell)
                           .real();
-          Ez_error(kk - 1, jj, ii) = Ez_interp - Ez_exact;
-          Ez_split_error(kk - 1, jj, ii) = Ez_split_interp - Ez_exact;
+          Ez_error(ii, jj, kk - 1) = Ez_interp - Ez_exact;
+          Ez_split_error(ii, jj, kk - 1) = Ez_split_interp - Ez_exact;
         }
       }
     }
@@ -291,7 +289,7 @@ TEST_CASE("E-field interpolation check") {
  *
  * Each component of the H-field will take the form
  * H_t(tt) = sin(2\pi tt) * exp(-tt^2)
- * The decision to make this function wholey imaginary is just to test whether
+ * The decision to make this function wholly imaginary is just to test whether
  * the imaginary part of the field is correctly picked up and worked with.
  *
  * We only test Fro-norm error metrics, since interpolation must occur along two