diff --git a/tdms/include/arrays/tensor3d.h b/tdms/include/arrays/tensor3d.h index e8db136a2..6c2861af5 100644 --- a/tdms/include/arrays/tensor3d.h +++ b/tdms/include/arrays/tensor3d.h @@ -22,18 +22,31 @@ * @tparam T Numerical datatype */ template -class Tensor3D : public std::vector { +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 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(){}; + Tensor3D() = default; Tensor3D(int n_layers, int n_cols, int n_rows) { allocate(n_layers, n_cols, n_rows); } @@ -43,20 +56,18 @@ class Tensor3D : public std::vector { /** @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)]; } /** @@ -64,7 +75,7 @@ class Tensor3D : public std::vector { * @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 { n_layers_ = n_layers; n_cols_ = n_cols; n_rows_ = n_rows; - this->resize(this->total_elements()); + data_.resize(total_elements()); } /** @@ -93,12 +104,12 @@ class Tensor3D : public std::vector { */ 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]; } } } @@ -106,7 +117,7 @@ class Tensor3D : public std::vector { 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]; } } } @@ -114,7 +125,7 @@ class Tensor3D : public std::vector { } /** @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 { */ 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); diff --git a/tdms/tests/unit/field_tests/test_Field_interpolation.cpp b/tdms/tests/unit/field_tests/test_Field_interpolation.cpp index d94f956fb..2af005c44 100644 --- a/tdms/tests/unit/field_tests/test_Field_interpolation.cpp +++ b/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 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 Ex_error(Nz, Ny, Nx - 1); + Tensor3D Ex_split_error(Nz, Ny, Nx - 1); + Tensor3D Ey_error(Nz, Ny - 1, Nx); + Tensor3D Ey_split_error(Nz, Ny - 1, Nx); + Tensor3D Ez_error(Nz - 1, Ny, Nx); + Tensor3D 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