diff --git a/tdms/include/arrays.h b/tdms/include/arrays.h index 22e1b7bdd..f428d49c4 100644 --- a/tdms/include/arrays.h +++ b/tdms/include/arrays.h @@ -11,6 +11,7 @@ #include +#include "arrays/tensor3d.h" #include "globals.h" #include "matlabio.h" #include "utils.h" @@ -354,127 +355,6 @@ class Pupil : public Matrix { ~Pupil(); }; -template -class Tensor3D { -protected: - int n_layers = 0; - int n_cols = 0; - int n_rows = 0; - T ***tensor = nullptr; - -public: - bool is_matlab_initialised = false; - - Tensor3D() = default; - - Tensor3D(T ***tensor, int n_layers, int n_cols, int n_rows) { - initialise(tensor, n_layers, n_cols, n_rows); - } - - void initialise(T ***_tensor, int _n_layers, int _n_cols, int _n_rows) { - tensor = _tensor; - n_layers = _n_layers; - n_cols = _n_cols; - n_rows = _n_rows; - } - - inline T **operator[](int value) const { return tensor[value]; }; - - bool has_elements() { return tensor != nullptr; }; - - void zero() { - for (int k = 0; k < n_layers; k++) - for (int j = 0; j < n_cols; j++) - for (int i = 0; i < n_rows; i++) { tensor[k][j][i] = 0; } - } - - void allocate(int nK, int nJ, int nI) { - n_layers = nK, n_cols = nJ, n_rows = nI; - tensor = (T ***) malloc(n_layers * sizeof(T **)); - - for (int k = 0; k < n_layers; k++) { - tensor[k] = (T **) malloc(n_cols * sizeof(T *)); - } - - for (int k = 0; k < n_layers; k++) { - for (int j = 0; j < n_cols; j++) { - tensor[k][j] = (T *) malloc(n_rows * sizeof(T)); - } - } - }; - - /** - * @brief Computes the Frobenius norm of the tensor - * - * fro_norm = \f$\sqrt{ \sum_{i=0}^{I_tot}\sum_{j=0}^{J_tot}\sum_{k=0}^{K_tot} - * |t[k][j][i]|^2 }\f$ - */ - double frobenius() { - T norm_val = 0; - for (int i1 = 0; i1 < n_layers; i1++) { - for (int i2 = 0; i2 < n_cols; i2++) { - for (int i3 = 0; i3 < n_rows; i3++) { - norm_val += abs(tensor[i1][i2][i3]) * abs(tensor[i1][i2][i3]); - } - } - } - return sqrt(norm_val); - } - - ~Tensor3D() { - if (tensor == nullptr) return; - if (is_matlab_initialised) { - free_cast_matlab_3D_array(tensor, n_layers); - } else { - for (int k = 0; k < n_layers; k++) { - for (int j = 0; j < n_cols; j++) { free(tensor[k][j]); } - free(tensor[k]); - } - free(tensor); - } - } -}; - -/** - * @brief Stores the fibre modes in the Fourier plane of the objective lens. - * - * The "Tilde" indicates that these quantities are in a Fourier plane relative - * to where the optical fibre is actually located, meaning that is has a Fourier - * relationship relative to the physical fibre mode(s). - */ -class DTilde { -protected: - int n_det_modes = 0;//< Number of modes specified - static void set_component(Tensor3D> &tensor, - const mxArray *ptr, const std::string &name, - int n_rows, int n_cols); - -public: - /** @brief Fetch the number of modes */ - inline int num_det_modes() const { return n_det_modes; }; - - /*! 3-dimensional vector of fibre modes indexed by (j, i, i_m). - * i and j index over the x and y plane respectively. - * i_m indexes the different modes specified in the input file.*/ - Tensor3D> x; - /*! @copydoc x */ - Tensor3D> y; - - void initialise(const mxArray *ptr, int n_rows, int n_cols); -}; - -class IncidentField { -protected: - void set_component(Tensor3D &component, const mxArray *ptr, - const std::string &name); - -public: - Tensor3D x; - Tensor3D y; - - explicit IncidentField(const mxArray *ptr); -}; - /** * List of field components as integers */ diff --git a/tdms/include/arrays/dtilde.h b/tdms/include/arrays/dtilde.h new file mode 100644 index 000000000..b25047ab3 --- /dev/null +++ b/tdms/include/arrays/dtilde.h @@ -0,0 +1,55 @@ +/** + * @file dtilde.h + * @brief Container class for the DTilde object, storing the fibre modes in the + * Fourier plane of the objective lens. + */ +#pragma once + +#include +#include +#include + +#include "arrays/tensor3d.h" +#include "matrix.h" + +/** + * @brief Stores the fibre modes in the Fourier plane of the objective lens. + * @details The "Tilde" indicates that these quantities are in a Fourier plane + * relative to where the optical fibre is actually located, meaning + * that it has a Fourier relationship relative to the physical fibre + * mode(s). + */ +class DTilde { +protected: + int n_det_modes = 0;//> &tensor, + const mxArray *ptr, const std::string &name, + int n_rows, int n_cols); + +public: + DTilde() = default; + + /** @brief (Fetch the) number of fibre modes specified */ + int num_det_modes() const { return n_det_modes; }; + + Tensor3D> x; + Tensor3D> y; + + /** + * @brief Initialise the fibre modes by reading from a MATLAB buffer. + * + * @param ptr Pointer to a MATLAB struct containing the data to read + * @param n_rows,n_cols Dimensions to assign to the component + */ + void initialise(const mxArray *ptr, int n_rows, int n_cols); +}; diff --git a/tdms/include/arrays/incident_field.h b/tdms/include/arrays/incident_field.h new file mode 100644 index 000000000..f207d37eb --- /dev/null +++ b/tdms/include/arrays/incident_field.h @@ -0,0 +1,22 @@ +/** + * @file incident_field.h + * @brief Declaration for the IncidentField object. + */ +#pragma once + +#include + +#include "arrays/tensor3d.h" +#include "matrix.h" + +class IncidentField { +protected: + void set_component(Tensor3D &component, const mxArray *ptr, + const std::string &name); + +public: + Tensor3D x; + Tensor3D y; + + explicit IncidentField(const mxArray *ptr); +}; diff --git a/tdms/include/arrays/tensor3d.h b/tdms/include/arrays/tensor3d.h new file mode 100644 index 000000000..5cc8f72a2 --- /dev/null +++ b/tdms/include/arrays/tensor3d.h @@ -0,0 +1,153 @@ +/** + * @file tensor3d.h + * @author William Graham + * @brief template class for storing three-dimensional data arrays. + */ +#pragma once + +#include + +#include "cell_coordinate.h" + +/** + * @brief Template class for storing three-dimensional data as a strided vector. + * @details Three dimensional data is stored as a strided vector, in row-major + * (C) format. The last listed dimension has the fastest varying index, and the + * first listed dimension has the slowest listed. Explicitly, + * + * this->(i, j, k) = this[i*n_cols*n_layers + j*n_layers + k]. + * + * This is consistent with the way hdf5 files store array-like data, see + * https://support.hdfgroup.org/HDF5/doc1.6/UG/12_Dataspaces.html. + * @tparam T Numerical datatype + */ +template +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() = default; + Tensor3D(int n_layers, int n_cols, int n_rows) { + allocate(n_layers, n_cols, n_rows); + } + Tensor3D(T ***buffer, int n_layers, int n_cols, int n_rows) { + initialise(buffer, n_layers, n_cols, n_rows); + } + + /** + * @brief Subscript operator for the Tensor, retrieves the (i,j,k)-th element. + * @note Does not 'call' the tensor, but rather accesses the data. + * @example ... + * @seealso Tensor3D::operator[](const ijk&) + */ + 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 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 data_[to_global_index(index_3d)]; + } + T operator[](const ijk &index_3d) const { + 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 total_elements() != 0; } + + /** + * @brief Allocate memory for this tensor given the dimensions passed. + * + * @param n_layers,n_cols,n_rows Number of layers, columns, and rows to + * assign. + */ + void allocate(int n_layers, int n_cols, int n_rows) { + n_layers_ = n_layers; + n_cols_ = n_cols; + n_rows_ = n_rows; + data_.resize(total_elements()); + } + + /** + * @brief Initialise this tensor from a 3D-buffer of matching size. + * @details Data values are copied so that membership over this->data() is + * preserved. + * @param buffer 3D buffer to read from + * @param n_layers,n_cols,n_rows "Shape" of the read buffer to assign to this + * tensor. + * @param buffer_leads_n_layers If true, the buffer to read from is + * assumed to have dimensions [n_layers][n_cols][n_rows]. If false, it is + * assumed to have the dimensions in reverse order, + * [n_rows][n_cols][n_layers]. + */ + void initialise(T ***buffer, int n_layers, int n_cols, int n_rows, + bool buffer_leads_n_layers = false) { + 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++) { + 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++) { + operator()(i, j, k) = buffer[i][j][k]; + } + } + } + } + } + + /** @brief Set all elements in the tensor to 0. */ + void zero() { std::fill(data_.begin(), data_.end(), 0); } + + /** + * @brief Computes the Frobenius norm of the tensor. + * @details The Frobenius norm is defined as: + * + * fro_norm = \f$\sqrt{ + * \sum_{i=0}^{n_rows}\sum_{j=0}^{n_cols}\sum_{k=0}^{n_layers} |t[k][j][i]|^2 + * }\f$. + * + * In practice, our data is flat and since addition is commutative, we can + * just iterate over the flattened data structure instead of using nested for + * loops. + */ + double frobenius() const { + T norm_val = 0; + 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/include/field.h b/tdms/include/field.h index f4cd2a69a..56d6fd6e5 100644 --- a/tdms/include/field.h +++ b/tdms/include/field.h @@ -61,9 +61,6 @@ class SplitFieldComponent : public Tensor3D { fftw_plan *plan_f = nullptr;// Forward fftw plan fftw_plan *plan_b = nullptr;// Backward fftw plan - double **operator[](int value) const { return tensor[value]; }; - double &operator[](ijk cell) const { return tensor[cell.k][cell.j][cell.i]; } - void initialise_from_matlab(double ***tensor, Dimensions &dims); /** diff --git a/tdms/include/simulation_manager/objects_from_infile.h b/tdms/include/simulation_manager/objects_from_infile.h index 22fde01aa..9d6fecfbe 100644 --- a/tdms/include/simulation_manager/objects_from_infile.h +++ b/tdms/include/simulation_manager/objects_from_infile.h @@ -11,6 +11,8 @@ #include "arrays.h" #include "arrays/cuboid.h" +#include "arrays/dtilde.h" +#include "arrays/incident_field.h" #include "cell_coordinate.h" #include "field.h" #include "fieldsample.h" diff --git a/tdms/src/arrays.cpp b/tdms/src/arrays.cpp index f94b5155d..e19782fd2 100644 --- a/tdms/src/arrays.cpp +++ b/tdms/src/arrays.cpp @@ -242,85 +242,6 @@ Pupil::~Pupil() { matrix = nullptr; } -void DTilde::set_component(Tensor3D> &tensor, - const mxArray *ptr, const string &name, int n_rows, - int n_cols) { - - auto element = ptr_to_nd_array_in(ptr, 3, name, "D_tilde"); - - auto dims = (int *) mxGetDimensions(element); - int n_det_modes = dims[0]; - - if (dims[1] != n_rows || dims[2] != n_cols) { - throw runtime_error("D_tilde.{x, y} has final dimensions " + - to_string(dims[1]) + "x" + to_string(dims[2]) + - " but it needed to be " + to_string(n_rows) + "x" + - to_string(n_cols)); - } - - auto p = (complex ***) malloc(sizeof(complex **) * n_cols); - for (int j = 0; j < n_cols; j++) { - p[j] = (complex **) malloc(sizeof(complex *) * n_rows); - for (int i = 0; i < n_rows; i++) { - p[j][i] = - (complex *) malloc(sizeof(complex) * n_det_modes); - } - } - - auto temp_re = - cast_matlab_3D_array(mxGetPr(element), dims[0], dims[1], dims[2]); - auto temp_im = - cast_matlab_3D_array(mxGetPi(element), dims[0], dims[1], dims[2]); - - for (int k = 0; k < n_det_modes; k++) - for (int j = 0; j < n_cols; j++) - for (int i = 0; i < n_rows; i++) { - p[j][i][k] = temp_re[j][i][k] + IMAGINARY_UNIT * temp_im[j][i][k]; - } - - free_cast_matlab_3D_array(temp_re, n_cols); - free_cast_matlab_3D_array(temp_im, n_cols); - tensor.initialise(p, n_cols, n_rows, n_det_modes); -} - -void DTilde::initialise(const mxArray *ptr, int n_rows, int n_cols) { - - if (mxIsEmpty(ptr)) { return; } - - assert_is_struct_with_n_fields(ptr, 2, "D_tilde"); - set_component(x, ptr, "Dx_tilde", n_rows, n_cols); - set_component(y, ptr, "Dy_tilde", n_rows, n_cols); - n_det_modes = - mxGetDimensions(ptr_to_nd_array_in(ptr, 3, "Dx_tilde", "D_tilde"))[0]; -} - -void IncidentField::set_component(Tensor3D &component, - const mxArray *ptr, const std::string &name) { - - if (mxIsEmpty(mxGetField(ptr, 0, name.c_str()))) { - spdlog::info("{} not present", name); - return; - } - - auto element = ptr_to_nd_array_in(ptr, 3, name, "tdfield"); - auto dims = mxGetDimensions(element); - int N = dims[0], M = dims[1], O = dims[2]; - component.initialise(cast_matlab_3D_array(mxGetPr(element), N, M, O), O, M, - N); - component.is_matlab_initialised = true; - - cerr << "Got tdfield, dims=(" + to_string(N) + "," + to_string(M) + "," + - to_string(O) + ")" - << endl; -} - -IncidentField::IncidentField(const mxArray *ptr) { - - assert_is_struct_with_n_fields(ptr, 2, "tdfield"); - set_component(x, ptr, "exi"); - set_component(y, ptr, "eyi"); -} - void FieldComponentsVector::initialise(const mxArray *ptr) { auto element = ptr_to_matrix_in(ptr, "components", "campssample"); diff --git a/tdms/src/arrays/dtilde.cpp b/tdms/src/arrays/dtilde.cpp new file mode 100644 index 000000000..88da92a99 --- /dev/null +++ b/tdms/src/arrays/dtilde.cpp @@ -0,0 +1,61 @@ +#include "arrays/dtilde.h" + +#include "globals.h" +#include "matlabio.h" +#include "matrix.h" + +using namespace std; +using tdms_math_constants::IMAGINARY_UNIT; + +void DTilde::set_component(Tensor3D> &tensor, + const mxArray *ptr, const string &name, int n_rows, + int n_cols) { + + auto element = ptr_to_nd_array_in(ptr, 3, name, "D_tilde"); + + auto dims = (int *) mxGetDimensions(element); + int n_det_modes = dims[0]; + + if (dims[1] != n_rows || dims[2] != n_cols) { + throw runtime_error("D_tilde.{x, y} has final dimensions " + + to_string(dims[1]) + "x" + to_string(dims[2]) + + " but it needed to be " + to_string(n_rows) + "x" + + to_string(n_cols)); + } + + complex ***p = + (complex ***) malloc(sizeof(complex **) * n_cols); + for (int j = 0; j < n_cols; j++) { + p[j] = (complex **) malloc(sizeof(complex *) * n_rows); + for (int i = 0; i < n_rows; i++) { + p[j][i] = + (complex *) malloc(sizeof(complex) * n_det_modes); + } + } + + auto temp_re = + cast_matlab_3D_array(mxGetPr(element), dims[0], dims[1], dims[2]); + auto temp_im = + cast_matlab_3D_array(mxGetPi(element), dims[0], dims[1], dims[2]); + + for (int k = 0; k < n_det_modes; k++) + for (int j = 0; j < n_cols; j++) + for (int i = 0; i < n_rows; i++) { + p[j][i][k] = temp_re[j][i][k] + IMAGINARY_UNIT * temp_im[j][i][k]; + } + + free_cast_matlab_3D_array(temp_re, n_cols); + free_cast_matlab_3D_array(temp_im, n_cols); + tensor.initialise(p, n_cols, n_rows, n_det_modes, true); +} + +void DTilde::initialise(const mxArray *ptr, int n_rows, int n_cols) { + + if (mxIsEmpty(ptr)) { return; } + + assert_is_struct_with_n_fields(ptr, 2, "D_tilde"); + set_component(x, ptr, "Dx_tilde", n_rows, n_cols); + set_component(y, ptr, "Dy_tilde", n_rows, n_cols); + n_det_modes = + mxGetDimensions(ptr_to_nd_array_in(ptr, 3, "Dx_tilde", "D_tilde"))[0]; +} diff --git a/tdms/src/arrays/indident_field.cpp b/tdms/src/arrays/indident_field.cpp new file mode 100644 index 000000000..7fffe7211 --- /dev/null +++ b/tdms/src/arrays/indident_field.cpp @@ -0,0 +1,31 @@ +#include "arrays/incident_field.h" + +#include + +#include "matlabio.h" + +using namespace std; + +void IncidentField::set_component(Tensor3D &component, + const mxArray *ptr, const std::string &name) { + + if (mxIsEmpty(mxGetField(ptr, 0, name.c_str()))) { + spdlog::info("{} not present", name); + return; + } + + auto element = ptr_to_nd_array_in(ptr, 3, name, "tdfield"); + auto dims = mxGetDimensions(element); + int N = dims[0], M = dims[1], O = dims[2]; + component.initialise(cast_matlab_3D_array(mxGetPr(element), N, M, O), O, M, N, + true); + + spdlog::info("Got tdfield, dims=({},{},{})", N, M, O); +} + +IncidentField::IncidentField(const mxArray *ptr) { + + assert_is_struct_with_n_fields(ptr, 2, "tdfield"); + set_component(x, ptr, "exi"); + set_component(y, ptr, "eyi"); +} diff --git a/tdms/src/fields/base.cpp b/tdms/src/fields/base.cpp index 32d8a8464..e6b3f2810 100644 --- a/tdms/src/fields/base.cpp +++ b/tdms/src/fields/base.cpp @@ -212,9 +212,9 @@ void Field::set_phasors(SplitField &F, int n, double omega, double dt, int Nt) { for (int j = jl; j <= ju; j++) for (int i = il; i <= iu; i++) { - x_m = F.xy[k][j][i] + F.xz[k][j][i]; - y_m = F.yx[k][j][i] + F.yz[k][j][i]; - z_m = F.zx[k][j][i] + F.zy[k][j][i]; + x_m = F.xy(i, j, k) + F.xz(i, j, k); + y_m = F.yx(i, j, k) + F.yz(i, j, k); + z_m = F.zx(i, j, k) + F.zy(i, j, k); int di = i - il; int dj = j - jl; diff --git a/tdms/src/fields/electric.cpp b/tdms/src/fields/electric.cpp index 49697113f..d88eaf7a3 100644 --- a/tdms/src/fields/electric.cpp +++ b/tdms/src/fields/electric.cpp @@ -113,8 +113,8 @@ double ElectricSplitField::interpolate_to_centre_of(AxialDirection d, for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { interp_data[ind] = - xy[k][j][i - scheme->number_of_datapoints_to_left + ind] + - xz[k][j][i - scheme->number_of_datapoints_to_left + ind]; + xy(i - scheme->number_of_datapoints_to_left + ind, j, k) + + xz(i - scheme->number_of_datapoints_to_left + ind, j, k); } // now run the interpolation scheme and place the result into the output return scheme->interpolate(interp_data); @@ -123,7 +123,7 @@ double ElectricSplitField::interpolate_to_centre_of(AxialDirection d, // if we are in a 2D simulation, we just return the field value at cell // (i, 0, k) since there is no y-dimension to interpolate in. if (tot.j <= 1) { - return yx[k][0][i] + yz[k][0][i]; + return yx(i, 0, k) + yz(i, 0, k); } else {// 3D simulation, interpolation is as normal scheme = &(best_scheme(tot.j, j, interpolation_method)); // now fill the interpolation data @@ -132,8 +132,8 @@ double ElectricSplitField::interpolate_to_centre_of(AxialDirection d, for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { interp_data[ind] = - yx[k][j - scheme->number_of_datapoints_to_left + ind][i] + - yz[k][j - scheme->number_of_datapoints_to_left + ind][i]; + yx(i, j - scheme->number_of_datapoints_to_left + ind, k) + + yz(i, j - scheme->number_of_datapoints_to_left + ind, k); } // now run the interpolation scheme and place the result into the output return scheme->interpolate(interp_data); @@ -146,8 +146,8 @@ double ElectricSplitField::interpolate_to_centre_of(AxialDirection d, for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { interp_data[ind] = - zx[k - scheme->number_of_datapoints_to_left + ind][j][i] + - zy[k - scheme->number_of_datapoints_to_left + ind][j][i]; + zx(i, j, k - scheme->number_of_datapoints_to_left + ind) + + zy(i, j, k - scheme->number_of_datapoints_to_left + ind); } // now run the interpolation scheme and place the result into the // output diff --git a/tdms/src/fields/magnetic.cpp b/tdms/src/fields/magnetic.cpp index e2685f5c7..a77a7abb7 100644 --- a/tdms/src/fields/magnetic.cpp +++ b/tdms/src/fields/magnetic.cpp @@ -341,8 +341,8 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, for (int ind = b_scheme->first_nonzero_coeff; ind <= b_scheme->last_nonzero_coeff; ind++) { data_for_first_scheme[ind] = - xy[k - b_scheme->number_of_datapoints_to_left + ind][j][i] + - xz[k - b_scheme->number_of_datapoints_to_left + ind][j][i]; + xy(i, j, k - b_scheme->number_of_datapoints_to_left + ind) + + xz(i, j, k - b_scheme->number_of_datapoints_to_left + ind); } // now run the interpolation scheme and place the result into the output @@ -368,7 +368,7 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int cell_k = k - c_scheme->number_of_datapoints_to_left + kk; // gather the data for interpolating in the z dimension data_for_first_scheme[kk] = - xy[cell_k][cell_j][i] + xz[cell_k][cell_j][i]; + xy(i, cell_j, cell_k) + xz(i, cell_j, cell_k); } // interpolate in z to obtain a value for the Hx field at position // (i, cell_j+Dy, k) place this into the appropriate index in the @@ -395,7 +395,7 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int cell_j = j - b_scheme->number_of_datapoints_to_left + jj; // gather the data for interpolating in the y dimension data_for_first_scheme[jj] = - xy[cell_k][cell_j][i] + xz[cell_k][cell_j][i]; + xy(i, cell_j, cell_k) + xz(i, cell_j, cell_k); } // interpolate in y to obtain a value for the Hx field at position // (i, j, cell_k+Dz) place this into the appropriate index in the @@ -432,7 +432,7 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int cell_k = k - b_scheme->number_of_datapoints_to_left + kk; // gather the data for interpolating in the z dimension data_for_first_scheme[kk] = - yx[cell_k][j][cell_i] + yz[cell_k][j][cell_i]; + yx(cell_i, j, cell_k) + yz(cell_i, j, cell_k); } // interpolate in z to obtain a value for the Hy field at position // (cell_i+Dx, j, k) place this into the appropriate index in the data @@ -458,7 +458,7 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int cell_i = i - c_scheme->number_of_datapoints_to_left + ii; // gather the data for interpolating in the x dimension data_for_first_scheme[ii] = - yx[cell_k][j][cell_i] + yz[cell_k][j][cell_i]; + yx(cell_i, j, cell_k) + yz(cell_i, j, cell_k); } // interpolate in x to obtain a value for the Hy field at position (i, // j, cell_k+Dz) place this into the appropriate index in the data @@ -481,8 +481,8 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, for (int ind = b_scheme->first_nonzero_coeff; ind <= b_scheme->last_nonzero_coeff; ind++) { data_for_first_scheme[ind] = - zx[k][j][i - b_scheme->number_of_datapoints_to_left + ind] + - zy[k][j][i - b_scheme->number_of_datapoints_to_left + ind]; + zx(i - b_scheme->number_of_datapoints_to_left + ind, j, k) + + zy(i - b_scheme->number_of_datapoints_to_left + ind, j, k); } // now run the interpolation scheme and place the result into the output @@ -508,7 +508,7 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int cell_j = j - c_scheme->number_of_datapoints_to_left + jj; // gather the data for interpolating in the y dimension data_for_first_scheme[jj] = - zx[k][cell_j][cell_i] + zy[k][cell_j][cell_i]; + zx(cell_i, cell_j, k) + zy(cell_i, cell_j, k); } // interpolate in y to obtain a value for the Hz field at position // (cell_i+Dx, j, k) place this into the appropriate index in the @@ -535,7 +535,7 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int cell_i = i - b_scheme->number_of_datapoints_to_left + ii; // gather the data for interpolating in the x dimension data_for_first_scheme[ii] = - zx[k][cell_j][cell_i] + zy[k][cell_j][cell_i]; + zx(cell_i, cell_j, k) + zy(cell_i, cell_j, k); } // interpolate in x to obtain a value for the Hz field at position // (i, j, cell_k+Dz) place this into the appropriate index in the diff --git a/tdms/src/fields/split.cpp b/tdms/src/fields/split.cpp index fd755e2fb..36e82e981 100644 --- a/tdms/src/fields/split.cpp +++ b/tdms/src/fields/split.cpp @@ -24,11 +24,7 @@ void SplitFieldComponent::initialise_fftw_plan(int n_threads, int size, void SplitFieldComponent::initialise_from_matlab(double ***tensor, Dimensions &dims) { - this->tensor = tensor; - this->n_layers = dims[2]; - this->n_cols = dims[1]; - this->n_rows = dims[0]; - this->is_matlab_initialised = true; + initialise(tensor, dims[2], dims[1], dims[0], true); } SplitFieldComponent::~SplitFieldComponent() { @@ -39,7 +35,6 @@ SplitFieldComponent::~SplitFieldComponent() { for (int i = 0; i < n_threads; i++) { fftw_destroy_plan(plan[i]); } free(plan); } - // superclass Tensor3D destructor removes tensor memory } void SplitField::allocate() { @@ -80,9 +75,9 @@ double SplitField::largest_field_value() { for (int k = 0; k < (tot.k + 1); k++) { for (int j = 0; j < (tot.j + 1); j++) { for (int i = 0; i < (tot.i + 1); i++) { - double x_field = fabs(xy[k][j][i] + xz[k][j][i]); - double y_field = fabs(yx[k][j][i] + yz[k][j][i]); - double z_field = fabs(zx[k][j][i] + zy[k][j][i]); + double x_field = fabs(xy(i, j, k) + xz(i, j, k)); + double y_field = fabs(yx(i, j, k) + yz(i, j, k)); + double z_field = fabs(zx(i, j, k) + zy(i, j, k)); if (largest_value < x_field) { largest_value = x_field; } if (largest_value < y_field) { largest_value = y_field; } if (largest_value < z_field) { largest_value = z_field; } diff --git a/tdms/src/fields/td_field_exporter_2d.cpp b/tdms/src/fields/td_field_exporter_2d.cpp index 4f29eeba8..35bbfeba6 100644 --- a/tdms/src/fields/td_field_exporter_2d.cpp +++ b/tdms/src/fields/td_field_exporter_2d.cpp @@ -35,7 +35,7 @@ void TDFieldExporter2D::export_field(SplitField &F, int stride, while (i < F.tot.i) { int k = 0; while (k < F.tot.k) { - array[k][i] = F.xy[k][0][i] + F.xz[k][0][i]; + array[k][i] = F.xy(i, 0, k) + F.xz(i, 0, k); k += stride; } i += stride; diff --git a/tdms/src/simulation_manager/execute_detector_subfunctions.cpp b/tdms/src/simulation_manager/execute_detector_subfunctions.cpp index dab04f6ea..851f6a3b3 100644 --- a/tdms/src/simulation_manager/execute_detector_subfunctions.cpp +++ b/tdms/src/simulation_manager/execute_detector_subfunctions.cpp @@ -41,11 +41,11 @@ void SimulationManager::compute_detector_functions(unsigned int tind, int m = j - inputs.params.pml.Dyl + (i - inputs.params.pml.Dxl) * (J_tot - inputs.params.pml.Dyu - inputs.params.pml.Dyl); - lv.Ex_t.v[m][0] = inputs.E_s.xy[inputs.params.k_det_obs][j][i] + - inputs.E_s.xz[inputs.params.k_det_obs][j][i]; + lv.Ex_t.v[m][0] = inputs.E_s.xy(i, j, inputs.params.k_det_obs) + + inputs.E_s.xz(i, j, inputs.params.k_det_obs); lv.Ex_t.v[m][1] = 0.; - lv.Ey_t.v[m][0] = inputs.E_s.yx[inputs.params.k_det_obs][j][i] + - inputs.E_s.yz[inputs.params.k_det_obs][j][i]; + lv.Ey_t.v[m][0] = inputs.E_s.yx(i, j, inputs.params.k_det_obs) + + inputs.E_s.yz(i, j, inputs.params.k_det_obs); lv.Ey_t.v[m][1] = 0.; } @@ -71,8 +71,8 @@ void SimulationManager::compute_detector_functions(unsigned int tind, j++) for (int i = 0; i < (I_tot - inputs.params.pml.Dxu - inputs.params.pml.Dxl); i++) { - lv.Ex_t.cm[j][i] *= inputs.pupil[j][i] * inputs.D_tilde.x[j][i][im]; - lv.Ey_t.cm[j][i] *= inputs.pupil[j][i] * inputs.D_tilde.y[j][i][im]; + lv.Ex_t.cm[j][i] *= inputs.pupil[j][i] * inputs.D_tilde.x(im, i, j); + lv.Ey_t.cm[j][i] *= inputs.pupil[j][i] * inputs.D_tilde.y(im, i, j); } /* Now iterate over each frequency we are extracting phasors at. diff --git a/tdms/src/simulation_manager/execute_simulation.cpp b/tdms/src/simulation_manager/execute_simulation.cpp index 35a29a0a8..3d4f29fd4 100644 --- a/tdms/src/simulation_manager/execute_simulation.cpp +++ b/tdms/src/simulation_manager/execute_simulation.cpp @@ -263,40 +263,40 @@ void SimulationManager::execute() { } - Enp1 = Ca * inputs.E_s.yx[k][j][i] + - Cb * (inputs.H_s.zx[k][j][i - 1] + - inputs.H_s.zy[k][j][i - 1] - - inputs.H_s.zx[k][j][i] - inputs.H_s.zy[k][j][i]); + Enp1 = Ca * inputs.E_s.yx(i, j, k) + + Cb * (inputs.H_s.zx(i - 1, j, k) + + inputs.H_s.zy(i - 1, j, k) - + inputs.H_s.zx(i, j, k) - inputs.H_s.zy(i, j, k)); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * loop_variables.E_nm1.yx[k][j][i] - + Enp1 += Cc * loop_variables.E_nm1.yx(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dx * ((1 + alpha_l) * - loop_variables.J_s.yx[k][j][i] + - beta_l * loop_variables.J_nm1.yx[k][j][i]); + loop_variables.J_s.yx(i, j, k) + + beta_l * loop_variables.J_nm1.yx(i, j, k)); if (loop_variables.is_conductive && rho) Enp1 += Cb * inputs.params.delta.dx * - loop_variables.J_c.yx[k][j][i]; + loop_variables.J_c.yx(i, j, k); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * loop_variables.J_s.yx[k][j][i] + - beta_l * loop_variables.J_nm1.yx[k][j][i] + + Jnp1 = alpha_l * loop_variables.J_s.yx(i, j, k) + + beta_l * loop_variables.J_nm1.yx(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - loop_variables.E_nm1.yx[k][j][i]); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.yx[k][j][i]; - loop_variables.E_nm1.yx[k][j][i] = inputs.E_s.yx[k][j][i]; - loop_variables.J_nm1.yx[k][j][i] = - loop_variables.J_s.yx[k][j][i]; - loop_variables.J_s.yx[k][j][i] = Jnp1; + (Enp1 - loop_variables.E_nm1.yx(i, j, k)); + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.yx(i, j, k); + loop_variables.E_nm1.yx(i, j, k) = inputs.E_s.yx(i, j, k); + loop_variables.J_nm1.yx(i, j, k) = + loop_variables.J_s.yx(i, j, k); + loop_variables.J_s.yx(i, j, k) = Jnp1; } if (loop_variables.is_conductive && rho) { - loop_variables.J_c.yx[k][j][i] -= - rho * (Enp1 + inputs.E_s.yx[k][j][i]); + loop_variables.J_c.yx(i, j, k) -= + rho * (Enp1 + inputs.E_s.yx(i, j, k)); } - inputs.E_s.yx[k][j][i] = Enp1; + inputs.E_s.yx(i, j, k) = Enp1; } // FDTD, E_s.yx } else { @@ -421,45 +421,45 @@ void SimulationManager::execute() { } - // Enp1 = Ca*E_s.yx[k][j][i]+Cb*(H_s.zx[k][j][i-1] + - // H_s.zy[k][j][i-1] - H_s.zx[k][j][i] - H_s.zy[k][j][i]); + // Enp1 = Ca*E_s.yx(i,j,k)+Cb*(H_s.zx[k][j][i-1] + + // H_s.zy[k][j][i-1] - H_s.zx(i,j,k) - H_s.zy(i,j,k)); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * loop_variables.E_nm1.yx[k][j][i] - + Enp1 += Cc * loop_variables.E_nm1.yx(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dx * ((1 + alpha_l) * - loop_variables.J_s.yx[k][j][i] + - beta_l * loop_variables.J_nm1.yx[k][j][i]); + loop_variables.J_s.yx(i, j, k) + + beta_l * loop_variables.J_nm1.yx(i, j, k)); if (loop_variables.is_conductive && rho) Enp1 += Cb * inputs.params.delta.dx * - loop_variables.J_c.yx[k][j][i]; + loop_variables.J_c.yx(i, j, k); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * loop_variables.J_s.yx[k][j][i] + - beta_l * loop_variables.J_nm1.yx[k][j][i] + + Jnp1 = alpha_l * loop_variables.J_s.yx(i, j, k) + + beta_l * loop_variables.J_nm1.yx(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - loop_variables.E_nm1.yx[k][j][i]); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.yx[k][j][i]; - loop_variables.E_nm1.yx[k][j][i] = inputs.E_s.yx[k][j][i]; - loop_variables.J_nm1.yx[k][j][i] = - loop_variables.J_s.yx[k][j][i]; - loop_variables.J_s.yx[k][j][i] = Jnp1; + (Enp1 - loop_variables.E_nm1.yx(i, j, k)); + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.yx(i, j, k); + loop_variables.E_nm1.yx(i, j, k) = inputs.E_s.yx(i, j, k); + loop_variables.J_nm1.yx(i, j, k) = + loop_variables.J_s.yx(i, j, k); + loop_variables.J_s.yx(i, j, k) = Jnp1; } if (loop_variables.is_conductive && rho) { - loop_variables.J_c.yx[k][j][i] -= - rho * (Enp1 + inputs.E_s.yx[k][j][i]); + loop_variables.J_c.yx(i, j, k) -= + rho * (Enp1 + inputs.E_s.yx(i, j, k)); } eh_vec[n][i][0] = - inputs.H_s.zx[k][j][i] + inputs.H_s.zy[k][j][i]; + inputs.H_s.zx(i, j, k) + inputs.H_s.zy(i, j, k); eh_vec[n][i][1] = 0.; PSTD.ca[n][i - 1] = Ca; PSTD.cb[n][i - 1] = Cb; } i = 0; - eh_vec[n][i][0] = inputs.H_s.zx[k][j][i] + inputs.H_s.zy[k][j][i]; + eh_vec[n][i][0] = inputs.H_s.zx(i, j, k) + inputs.H_s.zy(i, j, k); eh_vec[n][i][1] = 0.; first_derivative(eh_vec[n], eh_vec[n], PSTD.dk_ex, PSTD.N_ex, @@ -467,11 +467,11 @@ void SimulationManager::execute() { inputs.E_s.yx.plan_b[n]); for (i = 1; i < I_tot; i++) { - inputs.E_s.yx[k][j][i] = - PSTD.ca[n][i - 1] * inputs.E_s.yx[k][j][i] - + inputs.E_s.yx(i, j, k) = + PSTD.ca[n][i - 1] * inputs.E_s.yx(i, j, k) - PSTD.cb[n][i - 1] * eh_vec[n][i][0] / ((double) PSTD.N_ex); - // E_s.yx[k][j][i] = Enp1; + // E_s.yx(i,j,k) = Enp1; } } // PSTD, E_s.yx @@ -594,41 +594,41 @@ void SimulationManager::execute() { gamma_l = gamma_l / 2.; } } - Enp1 = Ca * inputs.E_s.yz[k][j][i] + - Cb * (inputs.H_s.xy[k][j][i] + inputs.H_s.xz[k][j][i] - - inputs.H_s.xy[k - 1][j][i] - - inputs.H_s.xz[k - 1][j][i]); + Enp1 = Ca * inputs.E_s.yz(i, j, k) + + Cb * (inputs.H_s.xy(i, j, k) + inputs.H_s.xz(i, j, k) - + inputs.H_s.xy(i, j, k - 1) - + inputs.H_s.xz(i, j, k - 1)); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * loop_variables.E_nm1.yz[k][j][i] - + Enp1 += Cc * loop_variables.E_nm1.yz(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dz * ((1 + alpha_l) * - loop_variables.J_s.yz[k][j][i] + - beta_l * loop_variables.J_nm1.yz[k][j][i]); + loop_variables.J_s.yz(i, j, k) + + beta_l * loop_variables.J_nm1.yz(i, j, k)); if (loop_variables.is_conductive && rho) Enp1 += Cb * inputs.params.delta.dz * - loop_variables.J_c.yz[k][j][i]; + loop_variables.J_c.yz(i, j, k); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * loop_variables.J_s.yz[k][j][i] + - beta_l * loop_variables.J_nm1.yz[k][j][i] + + Jnp1 = alpha_l * loop_variables.J_s.yz(i, j, k) + + beta_l * loop_variables.J_nm1.yz(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - loop_variables.E_nm1.yz[k][j][i]); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.yz[k][j][i]; - loop_variables.E_nm1.yz[k][j][i] = inputs.E_s.yz[k][j][i]; - loop_variables.J_nm1.yz[k][j][i] = - loop_variables.J_s.yz[k][j][i]; - loop_variables.J_s.yz[k][j][i] = Jnp1; + (Enp1 - loop_variables.E_nm1.yz(i, j, k)); + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.yz(i, j, k); + loop_variables.E_nm1.yz(i, j, k) = inputs.E_s.yz(i, j, k); + loop_variables.J_nm1.yz(i, j, k) = + loop_variables.J_s.yz(i, j, k); + loop_variables.J_s.yz(i, j, k) = Jnp1; } if (loop_variables.is_conductive && rho) { - loop_variables.J_c.yz[k][j][i] -= - rho * (Enp1 + inputs.E_s.yz[k][j][i]); + loop_variables.J_c.yz(i, j, k) -= + rho * (Enp1 + inputs.E_s.yz(i, j, k)); } - inputs.E_s.yz[k][j][i] = Enp1; + inputs.E_s.yz(i, j, k) = Enp1; } // FDTD, E_s.yz } else { @@ -746,46 +746,46 @@ void SimulationManager::execute() { gamma_l = gamma_l / 2.; } } - // Enp1 = Ca*E_s.yz[k][j][i]+Cb*(H_s.xy[k][j][i] + - // H_s.xz[k][j][i] - H_s.xy[k-1][j][i] - H_s.xz[k-1][j][i]); + // Enp1 = Ca*E_s.yz(i,j,k)+Cb*(H_s.xy(i,j,k) + + // H_s.xz(i, j, k) - H_s.xy[k-1][j][i] - H_s.xz[k-1][j][i]); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * loop_variables.E_nm1.yz[k][j][i] - + Enp1 += Cc * loop_variables.E_nm1.yz(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dz * ((1 + alpha_l) * - loop_variables.J_s.yz[k][j][i] + - beta_l * loop_variables.J_nm1.yz[k][j][i]); + loop_variables.J_s.yz(i, j, k) + + beta_l * loop_variables.J_nm1.yz(i, j, k)); if (loop_variables.is_conductive && rho) Enp1 += Cb * inputs.params.delta.dz * - loop_variables.J_c.yz[k][j][i]; + loop_variables.J_c.yz(i, j, k); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * loop_variables.J_s.yz[k][j][i] + - beta_l * loop_variables.J_nm1.yz[k][j][i] + + Jnp1 = alpha_l * loop_variables.J_s.yz(i, j, k) + + beta_l * loop_variables.J_nm1.yz(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - loop_variables.E_nm1.yz[k][j][i]); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.yz[k][j][i]; - loop_variables.E_nm1.yz[k][j][i] = inputs.E_s.yz[k][j][i]; - loop_variables.J_nm1.yz[k][j][i] = - loop_variables.J_s.yz[k][j][i]; - loop_variables.J_s.yz[k][j][i] = Jnp1; + (Enp1 - loop_variables.E_nm1.yz(i, j, k)); + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.yz(i, j, k); + loop_variables.E_nm1.yz(i, j, k) = inputs.E_s.yz(i, j, k); + loop_variables.J_nm1.yz(i, j, k) = + loop_variables.J_s.yz(i, j, k); + loop_variables.J_s.yz(i, j, k) = Jnp1; } if (loop_variables.is_conductive && rho) { - loop_variables.J_c.yz[k][j][i] -= - rho * (Enp1 + inputs.E_s.yz[k][j][i]); + loop_variables.J_c.yz(i, j, k) -= + rho * (Enp1 + inputs.E_s.yz(i, j, k)); } eh_vec[n][k][0] = - inputs.H_s.xy[k][j][i] + inputs.H_s.xz[k][j][i]; + inputs.H_s.xy(i, j, k) + inputs.H_s.xz(i, j, k); eh_vec[n][k][1] = 0.; PSTD.ca[n][k - 1] = Ca; PSTD.cb[n][k - 1] = Cb; } k = 0; - eh_vec[n][k][0] = inputs.H_s.xy[k][j][i] + inputs.H_s.xz[k][j][i]; + eh_vec[n][k][0] = inputs.H_s.xy(i, j, k) + inputs.H_s.xz(i, j, k); eh_vec[n][k][1] = 0.; first_derivative(eh_vec[n], eh_vec[n], PSTD.dk_ez, PSTD.N_ez, inputs.E_s.yz.plan_f[n], @@ -793,11 +793,11 @@ void SimulationManager::execute() { for (k = 1; k < K_tot; k++) { - inputs.E_s.yz[k][j][i] = - PSTD.ca[n][k - 1] * inputs.E_s.yz[k][j][i] + + inputs.E_s.yz(i, j, k) = + PSTD.ca[n][k - 1] * inputs.E_s.yz(i, j, k) + PSTD.cb[n][k - 1] * eh_vec[n][k][0] / ((double) PSTD.N_ez); - // E_s.yz[k][j][i] = Enp1; + // E_s.yz(i,j,k) = Enp1; } } // PSTD, E_s.yz @@ -922,40 +922,40 @@ void SimulationManager::execute() { gamma_l = gamma_l / 2.; } } - Enp1 = Ca * inputs.E_s.zx[k][j][i] + - Cb * (inputs.H_s.yx[k][j][i] + inputs.H_s.yz[k][j][i] - - inputs.H_s.yx[k][j][i - 1] - - inputs.H_s.yz[k][j][i - 1]); + Enp1 = Ca * inputs.E_s.zx(i, j, k) + + Cb * (inputs.H_s.yx(i, j, k) + inputs.H_s.yz(i, j, k) - + inputs.H_s.yx(i - 1, j, k) - + inputs.H_s.yz(i - 1, j, k)); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * loop_variables.E_nm1.zx[k][j][i] - + Enp1 += Cc * loop_variables.E_nm1.zx(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dx * ((1 + alpha_l) * - loop_variables.J_s.zx[k][j][i] + - beta_l * loop_variables.J_nm1.zx[k][j][i]); + loop_variables.J_s.zx(i, j, k) + + beta_l * loop_variables.J_nm1.zx(i, j, k)); if (loop_variables.is_conductive && rho) Enp1 += Cb * inputs.params.delta.dx * - loop_variables.J_c.zx[k][j][i]; + loop_variables.J_c.zx(i, j, k); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * loop_variables.J_s.zx[k][j][i] + - beta_l * loop_variables.J_nm1.zx[k][j][i] + + Jnp1 = alpha_l * loop_variables.J_s.zx(i, j, k) + + beta_l * loop_variables.J_nm1.zx(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - loop_variables.E_nm1.zx[k][j][i]); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zx[k][j][i]; - loop_variables.E_nm1.zx[k][j][i] = inputs.E_s.zx[k][j][i]; - loop_variables.J_nm1.zx[k][j][i] = - loop_variables.J_s.zx[k][j][i]; - loop_variables.J_s.zx[k][j][i] = Jnp1; + (Enp1 - loop_variables.E_nm1.zx(i, j, k)); + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zx(i, j, k); + loop_variables.E_nm1.zx(i, j, k) = inputs.E_s.zx(i, j, k); + loop_variables.J_nm1.zx(i, j, k) = + loop_variables.J_s.zx(i, j, k); + loop_variables.J_s.zx(i, j, k) = Jnp1; } if (loop_variables.is_conductive && rho) { - loop_variables.J_c.zx[k][j][i] -= - rho * (Enp1 + inputs.E_s.zx[k][j][i]); + loop_variables.J_c.zx(i, j, k) -= + rho * (Enp1 + inputs.E_s.zx(i, j, k)); } - inputs.E_s.zx[k][j][i] = Enp1; + inputs.E_s.zx(i, j, k) = Enp1; } // FDTD, E_s.zx } else { @@ -1072,45 +1072,45 @@ void SimulationManager::execute() { gamma_l = gamma_l / 2.; } } - // Enp1 = Ca*E_s.zx[k][j][i]+Cb*(H_s.yx[k][j][i] + - // H_s.yz[k][j][i] - H_s.yx[k][j][i-1] - H_s.yz[k][j][i-1]); + // Enp1 = Ca*E_s.zx(i,j,k)+Cb*(H_s.yx(i, j, k) + + // H_s.yz(i,j,k) - H_s.yx[k][j][i-1] - H_s.yz[k][j][i-1]); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * loop_variables.E_nm1.zx[k][j][i] - + Enp1 += Cc * loop_variables.E_nm1.zx(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dx * ((1 + alpha_l) * - loop_variables.J_s.zx[k][j][i] + - beta_l * loop_variables.J_nm1.zx[k][j][i]); + loop_variables.J_s.zx(i, j, k) + + beta_l * loop_variables.J_nm1.zx(i, j, k)); if (loop_variables.is_conductive && rho) Enp1 += Cb * inputs.params.delta.dx * - loop_variables.J_c.zx[k][j][i]; + loop_variables.J_c.zx(i, j, k); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * loop_variables.J_s.zx[k][j][i] + - beta_l * loop_variables.J_nm1.zx[k][j][i] + + Jnp1 = alpha_l * loop_variables.J_s.zx(i, j, k) + + beta_l * loop_variables.J_nm1.zx(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - loop_variables.E_nm1.zx[k][j][i]); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zx[k][j][i]; - loop_variables.E_nm1.zx[k][j][i] = inputs.E_s.zx[k][j][i]; - loop_variables.J_nm1.zx[k][j][i] = - loop_variables.J_s.zx[k][j][i]; - loop_variables.J_s.zx[k][j][i] = Jnp1; + (Enp1 - loop_variables.E_nm1.zx(i, j, k)); + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zx(i, j, k); + loop_variables.E_nm1.zx(i, j, k) = inputs.E_s.zx(i, j, k); + loop_variables.J_nm1.zx(i, j, k) = + loop_variables.J_s.zx(i, j, k); + loop_variables.J_s.zx(i, j, k) = Jnp1; } if (loop_variables.is_conductive && rho) { - loop_variables.J_c.zx[k][j][i] -= - rho * (Enp1 + inputs.E_s.zx[k][j][i]); + loop_variables.J_c.zx(i, j, k) -= + rho * (Enp1 + inputs.E_s.zx(i, j, k)); } eh_vec[n][i][0] = - inputs.H_s.yx[k][j][i] + inputs.H_s.yz[k][j][i]; + inputs.H_s.yx(i, j, k) + inputs.H_s.yz(i, j, k); eh_vec[n][i][1] = 0.; PSTD.ca[n][i - 1] = Ca; PSTD.cb[n][i - 1] = Cb; } i = 0; - eh_vec[n][i][0] = inputs.H_s.yx[k][j][i] + inputs.H_s.yz[k][j][i]; + eh_vec[n][i][0] = inputs.H_s.yx(i, j, k) + inputs.H_s.yz(i, j, k); eh_vec[n][i][1] = 0.; first_derivative(eh_vec[n], eh_vec[n], PSTD.dk_ex, PSTD.N_ex, @@ -1118,11 +1118,11 @@ void SimulationManager::execute() { inputs.E_s.zx.plan_b[n]); for (i = 1; i < I_tot; i++) { - inputs.E_s.zx[k][j][i] = - PSTD.ca[n][i - 1] * inputs.E_s.zx[k][j][i] + + inputs.E_s.zx(i, j, k) = + PSTD.ca[n][i - 1] * inputs.E_s.zx(i, j, k) + PSTD.cb[n][i - 1] * eh_vec[n][i][0] / ((double) PSTD.N_ex); - // E_s.zx[k][j][i] = Enp1; + // E_s.zx(i,j,k) = Enp1; } } // PSTD, E_s.zx @@ -1199,39 +1199,39 @@ void SimulationManager::execute() { } } - Enp1 = Ca * inputs.E_s.zx[k][j][i] + - Cb * (inputs.H_s.yx[k][j][i] + inputs.H_s.yz[k][j][i] - - inputs.H_s.yx[k][j][i - 1] - - inputs.H_s.yz[k][j][i - 1]); + Enp1 = Ca * inputs.E_s.zx(i, j, k) + + Cb * (inputs.H_s.yx(i, j, k) + inputs.H_s.yz(i, j, k) - + inputs.H_s.yx(i - 1, j, k) - + inputs.H_s.yz(i - 1, j, k)); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * loop_variables.E_nm1.zx[k][j][i] - + Enp1 += Cc * loop_variables.E_nm1.zx(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dx * ((1 + alpha_l) * - loop_variables.J_s.zx[k][j][i] + - beta_l * loop_variables.J_nm1.zx[k][j][i]); + loop_variables.J_s.zx(i, j, k) + + beta_l * loop_variables.J_nm1.zx(i, j, k)); if (loop_variables.is_conductive && rho) Enp1 += Cb * inputs.params.delta.dx * - loop_variables.J_c.zx[k][j][i]; + loop_variables.J_c.zx(i, j, k); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * loop_variables.J_s.zx[k][j][i] + - beta_l * loop_variables.J_nm1.zx[k][j][i] + + Jnp1 = alpha_l * loop_variables.J_s.zx(i, j, k) + + beta_l * loop_variables.J_nm1.zx(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - loop_variables.E_nm1.zx[k][j][i]); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zx[k][j][i]; - loop_variables.E_nm1.zx[k][j][i] = inputs.E_s.zx[k][j][i]; - loop_variables.J_nm1.zx[k][j][i] = - loop_variables.J_s.zx[k][j][i]; - loop_variables.J_s.zx[k][j][i] = Jnp1; + (Enp1 - loop_variables.E_nm1.zx(i, j, k)); + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zx(i, j, k); + loop_variables.E_nm1.zx(i, j, k) = inputs.E_s.zx(i, j, k); + loop_variables.J_nm1.zx(i, j, k) = + loop_variables.J_s.zx(i, j, k); + loop_variables.J_s.zx(i, j, k) = Jnp1; } if (loop_variables.is_conductive && rho) { - loop_variables.J_c.zx[k][j][i] -= - rho * (Enp1 + inputs.E_s.zx[k][j][i]); + loop_variables.J_c.zx(i, j, k) -= + rho * (Enp1 + inputs.E_s.zx(i, j, k)); } - inputs.E_s.zx[k][j][i] = Enp1; + inputs.E_s.zx(i, j, k) = Enp1; } } if (inputs.params.dimension == THREE || @@ -1353,41 +1353,41 @@ void SimulationManager::execute() { } - Enp1 = Ca * inputs.E_s.zy[k][j][i] + - Cb * (inputs.H_s.xy[k][j - 1][i] + - inputs.H_s.xz[k][j - 1][i] - - inputs.H_s.xy[k][j][i] - inputs.H_s.xz[k][j][i]); + Enp1 = Ca * inputs.E_s.zy(i, j, k) + + Cb * (inputs.H_s.xy(i, j - 1, k) + + inputs.H_s.xz(i, j - 1, k) - + inputs.H_s.xy(i, j, k) - inputs.H_s.xz(i, j, k)); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * loop_variables.E_nm1.zy[k][j][i] - + Enp1 += Cc * loop_variables.E_nm1.zy(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dy * ((1 + alpha_l) * - loop_variables.J_s.zy[k][j][i] + - beta_l * loop_variables.J_nm1.zy[k][j][i]); + loop_variables.J_s.zy(i, j, k) + + beta_l * loop_variables.J_nm1.zy(i, j, k)); if (loop_variables.is_conductive && rho) Enp1 += Cb * inputs.params.delta.dy * - loop_variables.J_c.zy[k][j][i]; + loop_variables.J_c.zy(i, j, k); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * loop_variables.J_s.zy[k][j][i] + - beta_l * loop_variables.J_nm1.zy[k][j][i] + + Jnp1 = alpha_l * loop_variables.J_s.zy(i, j, k) + + beta_l * loop_variables.J_nm1.zy(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - loop_variables.E_nm1.zy[k][j][i]); + (Enp1 - loop_variables.E_nm1.zy(i, j, k)); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zy[k][j][i]; - loop_variables.E_nm1.zy[k][j][i] = inputs.E_s.zy[k][j][i]; - loop_variables.J_nm1.zy[k][j][i] = - loop_variables.J_s.zy[k][j][i]; - loop_variables.J_s.zy[k][j][i] = Jnp1; + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zy(i, j, k); + loop_variables.E_nm1.zy(i, j, k) = inputs.E_s.zy(i, j, k); + loop_variables.J_nm1.zy(i, j, k) = + loop_variables.J_s.zy(i, j, k); + loop_variables.J_s.zy(i, j, k) = Jnp1; } if (loop_variables.is_conductive && rho) { - loop_variables.J_c.zy[k][j][i] -= - rho * (Enp1 + inputs.E_s.zy[k][j][i]); + loop_variables.J_c.zy(i, j, k) -= + rho * (Enp1 + inputs.E_s.zy(i, j, k)); } - inputs.E_s.zy[k][j][i] = Enp1; + inputs.E_s.zy(i, j, k) = Enp1; } // FDTD, E_s.zy } else { @@ -1506,41 +1506,41 @@ void SimulationManager::execute() { } - // Enp1 = Ca*E_s.zy[k][j][i]+Cb*(H_s.xy[k][j-1][i] + - // H_s.xz[k][j-1][i] - H_s.xy[k][j][i] - H_s.xz[k][j][i]); + // Enp1 = Ca*E_s.zy(i,j,k)+Cb*(H_s.xy[k][j-1][i] + + // H_s.xz[k][j-1][i] - H_s.xy(i,j,k) - H_s.xz(i, j, k)); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * loop_variables.E_nm1.zy[k][j][i] - + Enp1 += Cc * loop_variables.E_nm1.zy(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dy * ((1 + alpha_l) * - loop_variables.J_s.zy[k][j][i] + - beta_l * loop_variables.J_nm1.zy[k][j][i]); + loop_variables.J_s.zy(i, j, k) + + beta_l * loop_variables.J_nm1.zy(i, j, k)); if (loop_variables.is_conductive && rho) Enp1 += Cb * inputs.params.delta.dy * - loop_variables.J_c.zy[k][j][i]; + loop_variables.J_c.zy(i, j, k); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * loop_variables.J_s.zy[k][j][i] + - beta_l * loop_variables.J_nm1.zy[k][j][i] + + Jnp1 = alpha_l * loop_variables.J_s.zy(i, j, k) + + beta_l * loop_variables.J_nm1.zy(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - loop_variables.E_nm1.zy[k][j][i]); + (Enp1 - loop_variables.E_nm1.zy(i, j, k)); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zy[k][j][i]; - loop_variables.E_nm1.zy[k][j][i] = inputs.E_s.zy[k][j][i]; - loop_variables.J_nm1.zy[k][j][i] = - loop_variables.J_s.zy[k][j][i]; - loop_variables.J_s.zy[k][j][i] = Jnp1; + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zy(i, j, k); + loop_variables.E_nm1.zy(i, j, k) = inputs.E_s.zy(i, j, k); + loop_variables.J_nm1.zy(i, j, k) = + loop_variables.J_s.zy(i, j, k); + loop_variables.J_s.zy(i, j, k) = Jnp1; } if (loop_variables.is_conductive && rho) { - loop_variables.J_c.zy[k][j][i] -= - rho * (Enp1 + inputs.E_s.zy[k][j][i]); + loop_variables.J_c.zy(i, j, k) -= + rho * (Enp1 + inputs.E_s.zy(i, j, k)); } eh_vec[n][j][0] = - inputs.H_s.xy[k][j][i] + inputs.H_s.xz[k][j][i]; + inputs.H_s.xy(i, j, k) + inputs.H_s.xz(i, j, k); eh_vec[n][j][1] = 0.; PSTD.ca[n][j - 1] = Ca; PSTD.cb[n][j - 1] = Cb; @@ -1548,18 +1548,18 @@ void SimulationManager::execute() { if (J_tot > 1) { j = 0; eh_vec[n][j][0] = - inputs.H_s.xy[k][j][i] + inputs.H_s.xz[k][j][i]; + inputs.H_s.xy(i, j, k) + inputs.H_s.xz(i, j, k); eh_vec[n][j][1] = 0.; first_derivative(eh_vec[n], eh_vec[n], PSTD.dk_ey, PSTD.N_ey, inputs.E_s.zy.plan_f[n], inputs.E_s.zy.plan_b[n]); } for (j = 1; j < J_tot; j++) { - inputs.E_s.zy[k][j][i] = - PSTD.ca[n][j - 1] * inputs.E_s.zy[k][j][i] - + inputs.E_s.zy(i, j, k) = + PSTD.ca[n][j - 1] * inputs.E_s.zy(i, j, k) - PSTD.cb[n][j - 1] * eh_vec[n][j][0] / ((double) PSTD.N_ey); - // E_s.zy[k][j][i] = Enp1; + // E_s.zy(i,j,k) = Enp1; } } // PSTD, E_s.zy @@ -1635,40 +1635,40 @@ void SimulationManager::execute() { } - Enp1 = Ca * inputs.E_s.zy[k][j][i] + - Cb * (inputs.H_s.xy[k][j - 1][i] + - inputs.H_s.xz[k][j - 1][i] - inputs.H_s.xy[k][j][i] - - inputs.H_s.xz[k][j][i]); + Enp1 = Ca * inputs.E_s.zy(i, j, k) + + Cb * (inputs.H_s.xy(i, j - 1, k) + + inputs.H_s.xz(i, j - 1, k) - inputs.H_s.xy(i, j, k) - + inputs.H_s.xz(i, j, k)); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * loop_variables.E_nm1.zy[k][j][i] - + Enp1 += Cc * loop_variables.E_nm1.zy(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dy * ((1 + alpha_l) * - loop_variables.J_s.zy[k][j][i] + - beta_l * loop_variables.J_nm1.zy[k][j][i]); + loop_variables.J_s.zy(i, j, k) + + beta_l * loop_variables.J_nm1.zy(i, j, k)); if (loop_variables.is_conductive && rho) Enp1 += Cb * inputs.params.delta.dy * - loop_variables.J_c.zy[k][j][i]; + loop_variables.J_c.zy(i, j, k); if ((loop_variables.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * loop_variables.J_s.zy[k][j][i] + - beta_l * loop_variables.J_nm1.zy[k][j][i] + + Jnp1 = alpha_l * loop_variables.J_s.zy(i, j, k) + + beta_l * loop_variables.J_nm1.zy(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - loop_variables.E_nm1.zy[k][j][i]); + (Enp1 - loop_variables.E_nm1.zy(i, j, k)); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zy[k][j][i]; - loop_variables.E_nm1.zy[k][j][i] = inputs.E_s.zy[k][j][i]; - loop_variables.J_nm1.zy[k][j][i] = - loop_variables.J_s.zy[k][j][i]; - loop_variables.J_s.zy[k][j][i] = Jnp1; + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.zy(i, j, k); + loop_variables.E_nm1.zy(i, j, k) = inputs.E_s.zy(i, j, k); + loop_variables.J_nm1.zy(i, j, k) = + loop_variables.J_s.zy(i, j, k); + loop_variables.J_s.zy(i, j, k) = Jnp1; } if (loop_variables.is_conductive && rho) { - loop_variables.J_c.zy[k][j][i] -= - rho * (Enp1 + inputs.E_s.zy[k][j][i]); + loop_variables.J_c.zy(i, j, k) -= + rho * (Enp1 + inputs.E_s.zy(i, j, k)); } - inputs.E_s.zy[k][j][i] = Enp1; + inputs.E_s.zy(i, j, k) = Enp1; } } }// end of parallel section @@ -1752,21 +1752,21 @@ void SimulationManager::execute() { } if (!inputs.materials[k][j][i]) - inputs.H_s.xz[k][j][i] = - inputs.D.a.z[k_loc] * inputs.H_s.xz[k][j][i] + - inputs.D.b.z[k_loc] * (inputs.E_s.yx[k + 1][j][i] + - inputs.E_s.yz[k + 1][j][i] - - inputs.E_s.yx[k][j][i] - - inputs.E_s.yz[k][j][i]); + inputs.H_s.xz(i, j, k) = + inputs.D.a.z[k_loc] * inputs.H_s.xz(i, j, k) + + inputs.D.b.z[k_loc] * (inputs.E_s.yx(i, j, k + 1) + + inputs.E_s.yz(i, j, k + 1) - + inputs.E_s.yx(i, j, k) - + inputs.E_s.yz(i, j, k)); else - inputs.H_s.xz[k][j][i] = + inputs.H_s.xz(i, j, k) = inputs.Dmaterial.a.z[inputs.materials[k][j][i] - 1] * - inputs.H_s.xz[k][j][i] + + inputs.H_s.xz(i, j, k) + inputs.Dmaterial.b.z[inputs.materials[k][j][i] - 1] * - (inputs.E_s.yx[k + 1][j][i] + - inputs.E_s.yz[k + 1][j][i] - - inputs.E_s.yx[k][j][i] - - inputs.E_s.yz[k][j][i]); + (inputs.E_s.yx(i, j, k + 1) + + inputs.E_s.yz(i, j, k + 1) - + inputs.E_s.yx(i, j, k) - + inputs.E_s.yz(i, j, k)); } // FDTD, H_s.xz } else { @@ -1797,25 +1797,27 @@ void SimulationManager::execute() { if (!inputs.materials[k][j][i]) { PSTD.ca[n][k] = inputs.D.a.z[k_loc]; PSTD.cb[n][k] = inputs.D.b.z[k_loc]; - // H_s.xz[k][j][i] = - // D.a.z[k_loc]*H_s.xz[k][j][i]+D.b.z[k_loc]*(E_s.yx[k+1][j][i] - // + E_s.yz[k+1][j][i] - E_s.yx[k][j][i] - E_s.yz[k][j][i]); + // H_s.xz(i, j, k) = + // D.a.z[k_loc]*H_s.xz(i, j, + // k)+D.b.z[k_loc]*(E_s.yx[k+1][j][i] + // + E_s.yz[k+1][j][i] - E_s.yx(i,j,k) - E_s.yz(i,j,k)); } else { PSTD.ca[n][k] = inputs.Dmaterial.a.z[inputs.materials[k][j][i] - 1]; PSTD.cb[n][k] = inputs.Dmaterial.b.z[inputs.materials[k][j][i] - 1]; - // H_s.xz[k][j][i] = - // Dmaterial.Da.z[materials[k][j][i]-1]*H_s.xz[k][j][i]+Dmaterial.Db.z[materials[k][j][i]-1]*(E_s.yx[k+1][j][i] - // + E_s.yz[k+1][j][i] - E_s.yx[k][j][i] - E_s.yz[k][j][i]); + // H_s.xz(i, j, k) = + // Dmaterial.Da.z[materials[k][j][i]-1]*H_s.xz(i, j, + // k)+Dmaterial.Db.z[materials[k][j][i]-1]*(E_s.yx[k+1][j][i] + // + E_s.yz[k+1][j][i] - E_s.yx(i,j,k) - E_s.yz(i,j,k)); } eh_vec[n][k][0] = - inputs.E_s.yx[k][j][i] + inputs.E_s.yz[k][j][i]; + inputs.E_s.yx(i, j, k) + inputs.E_s.yz(i, j, k); eh_vec[n][k][1] = 0.; } k = K_tot; - eh_vec[n][k][0] = inputs.E_s.yx[k][j][i] + inputs.E_s.yz[k][j][i]; + eh_vec[n][k][0] = inputs.E_s.yx(i, j, k) + inputs.E_s.yz(i, j, k); eh_vec[n][k][1] = 0.; first_derivative(eh_vec[n], eh_vec[n], PSTD.dk_hz, PSTD.N_hz, @@ -1823,8 +1825,8 @@ void SimulationManager::execute() { inputs.H_s.xz.plan_b[n]); for (k = 0; k < K_tot; k++) { - inputs.H_s.xz[k][j][i] = - PSTD.ca[n][k] * inputs.H_s.xz[k][j][i] + + inputs.H_s.xz(i, j, k) = + PSTD.ca[n][k] * inputs.H_s.xz(i, j, k) + PSTD.cb[n][k] * eh_vec[n][k][0] / ((double) PSTD.N_hz); } } @@ -1862,22 +1864,22 @@ void SimulationManager::execute() { else array_ind = (J_tot + 1) * k_loc + j; if (!inputs.materials[k][j][i]) - inputs.H_s.xy[k][j][i] = - inputs.D.a.y[array_ind] * inputs.H_s.xy[k][j][i] + + inputs.H_s.xy(i, j, k) = + inputs.D.a.y[array_ind] * inputs.H_s.xy(i, j, k) + inputs.D.b.y[array_ind] * - (inputs.E_s.zy[k][j][i] + - inputs.E_s.zx[k][j][i] - - inputs.E_s.zy[k][j + 1][i] - - inputs.E_s.zx[k][j + 1][i]); + (inputs.E_s.zy(i, j, k) + + inputs.E_s.zx(i, j, k) - + inputs.E_s.zy(i, j + 1, k) - + inputs.E_s.zx(i, j + 1, k)); else - inputs.H_s.xy[k][j][i] = + inputs.H_s.xy(i, j, k) = inputs.Dmaterial.a.y[inputs.materials[k][j][i] - 1] * - inputs.H_s.xy[k][j][i] + + inputs.H_s.xy(i, j, k) + inputs.Dmaterial.b.y[inputs.materials[k][j][i] - 1] * - (inputs.E_s.zy[k][j][i] + - inputs.E_s.zx[k][j][i] - - inputs.E_s.zy[k][j + 1][i] - - inputs.E_s.zx[k][j + 1][i]); + (inputs.E_s.zy(i, j, k) + + inputs.E_s.zx(i, j, k) - + inputs.E_s.zy(i, j + 1, k) - + inputs.E_s.zx(i, j + 1, k)); } // FDTD, H_s.xy } else { @@ -1910,25 +1912,25 @@ void SimulationManager::execute() { if (!inputs.materials[k][j][i]) { PSTD.ca[n][j] = inputs.D.a.y[array_ind]; PSTD.cb[n][j] = inputs.D.b.y[array_ind]; - // H_s.xy[k][j][i] = - // D.a.y[array_ind]*H_s.xy[k][j][i]+D.b.y[array_ind]*(E_s.zy[k][j][i] - // + E_s.zx[k][j][i] - E_s.zy[k][j+1][i] - E_s.zx[k][j+1][i]); + // H_s.xy(i,j,k) = + // D.a.y[array_ind]*H_s.xy(i,j,k)+D.b.y[array_ind]*(E_s.zy(i,j,k) + // + E_s.zx(i,j,k) - E_s.zy[k][j+1][i] - E_s.zx[k][j+1][i]); } else { PSTD.ca[n][j] = inputs.Dmaterial.a.y[inputs.materials[k][j][i] - 1]; PSTD.cb[n][j] = inputs.Dmaterial.b.y[inputs.materials[k][j][i] - 1]; - // H_s.xy[k][j][i] = - // Dmaterial.Da.y[materials[k][j][i]-1]*H_s.xy[k][j][i]+Dmaterial.Db.y[materials[k][j][i]-1]*(E_s.zy[k][j][i] - //+ E_s.zx[k][j][i] - E_s.zy[k][j+1][i] - E_s.zx[k][j+1][i]); + // H_s.xy(i,j,k) = + // Dmaterial.Da.y[materials[k][j][i]-1]*H_s.xy(i,j,k)+Dmaterial.Db.y[materials[k][j][i]-1]*(E_s.zy(i,j,k) + //+ E_s.zx(i,j,k) - E_s.zy[k][j+1][i] - E_s.zx[k][j+1][i]); } eh_vec[n][j][0] = - inputs.E_s.zy[k][j][i] + inputs.E_s.zx[k][j][i]; + inputs.E_s.zy(i, j, k) + inputs.E_s.zx(i, j, k); eh_vec[n][j][1] = 0.; } j = J_tot; - eh_vec[n][j][0] = inputs.E_s.zy[k][j][i] + inputs.E_s.zx[k][j][i]; + eh_vec[n][j][0] = inputs.E_s.zy(i, j, k) + inputs.E_s.zx(i, j, k); eh_vec[n][j][1] = 0.; first_derivative(eh_vec[n], eh_vec[n], PSTD.dk_hy, PSTD.N_hy, @@ -1936,8 +1938,8 @@ void SimulationManager::execute() { inputs.H_s.xy.plan_b[n]); for (j = 0; j < J_tot; j++) { - inputs.H_s.xy[k][j][i] = - PSTD.ca[n][j] * inputs.H_s.xy[k][j][i] - + inputs.H_s.xy(i, j, k) = + PSTD.ca[n][j] * inputs.H_s.xy(i, j, k) - PSTD.cb[n][j] * eh_vec[n][j][0] / ((double) PSTD.N_hy); } } @@ -1974,22 +1976,22 @@ void SimulationManager::execute() { else array_ind = (I_tot + 1) * k_loc + i; if (!inputs.materials[k][j][i]) - inputs.H_s.yx[k][j][i] = - inputs.D.a.x[array_ind] * inputs.H_s.yx[k][j][i] + + inputs.H_s.yx(i, j, k) = + inputs.D.a.x[array_ind] * inputs.H_s.yx(i, j, k) + inputs.D.b.x[array_ind] * - (inputs.E_s.zx[k][j][i + 1] + - inputs.E_s.zy[k][j][i + 1] - - inputs.E_s.zx[k][j][i] - - inputs.E_s.zy[k][j][i]); + (inputs.E_s.zx(i + 1, j, k) + + inputs.E_s.zy(i + 1, j, k) - + inputs.E_s.zx(i, j, k) - + inputs.E_s.zy(i, j, k)); else { - inputs.H_s.yx[k][j][i] = + inputs.H_s.yx(i, j, k) = inputs.Dmaterial.a.x[inputs.materials[k][j][i] - 1] * - inputs.H_s.yx[k][j][i] + + inputs.H_s.yx(i, j, k) + inputs.Dmaterial.b.x[inputs.materials[k][j][i] - 1] * - (inputs.E_s.zx[k][j][i + 1] + - inputs.E_s.zy[k][j][i + 1] - - inputs.E_s.zx[k][j][i] - - inputs.E_s.zy[k][j][i]); + (inputs.E_s.zx(i + 1, j, k) + + inputs.E_s.zy(i + 1, j, k) - + inputs.E_s.zx(i, j, k) - + inputs.E_s.zy(i, j, k)); } } // FDTD, H_s.yx @@ -2023,25 +2025,27 @@ void SimulationManager::execute() { if (!inputs.materials[k][j][i]) { PSTD.ca[n][i] = inputs.D.a.x[array_ind]; PSTD.cb[n][i] = inputs.D.b.x[array_ind]; - // H_s.yx[k][j][i] = - // D.a.x[array_ind]*H_s.yx[k][j][i]+D.b.x[array_ind]*(E_s.zx[k][j][i+1] - // + E_s.zy[k][j][i+1] - E_s.zx[k][j][i] - E_s.zy[k][j][i]); + // H_s.yx(i, j, k) = + // D.a.x[array_ind]*H_s.yx(i, j, + // k)+D.b.x[array_ind]*(E_s.zx[k][j][i+1] + // + E_s.zy[k][j][i+1] - E_s.zx(i,j,k) - E_s.zy(i,j,k)); } else { PSTD.ca[n][i] = inputs.Dmaterial.a.x[inputs.materials[k][j][i] - 1]; PSTD.cb[n][i] = inputs.Dmaterial.b.x[inputs.materials[k][j][i] - 1]; - // H_s.yx[k][j][i] = - // Dmaterial.Da.x[materials[k][j][i]-1]*H_s.yx[k][j][i]+Dmaterial.Db.x[materials[k][j][i]-1]*(E_s.zx[k][j][i+1] - //+ E_s.zy[k][j][i+1] - E_s.zx[k][j][i] - E_s.zy[k][j][i]); + // H_s.yx(i, j, k) = + // Dmaterial.Da.x[materials[k][j][i]-1]*H_s.yx(i, j, + // k)+Dmaterial.Db.x[materials[k][j][i]-1]*(E_s.zx[k][j][i+1] + //+ E_s.zy[k][j][i+1] - E_s.zx(i,j,k) - E_s.zy(i,j,k)); } eh_vec[n][i][0] = - inputs.E_s.zx[k][j][i] + inputs.E_s.zy[k][j][i]; + inputs.E_s.zx(i, j, k) + inputs.E_s.zy(i, j, k); eh_vec[n][i][1] = 0.; } i = I_tot; - eh_vec[n][i][0] = inputs.E_s.zx[k][j][i] + inputs.E_s.zy[k][j][i]; + eh_vec[n][i][0] = inputs.E_s.zx(i, j, k) + inputs.E_s.zy(i, j, k); eh_vec[n][i][1] = 0.; first_derivative(eh_vec[n], eh_vec[n], PSTD.dk_hx, PSTD.N_hx, @@ -2049,8 +2053,8 @@ void SimulationManager::execute() { inputs.H_s.yx.plan_b[n]); for (i = 0; i < I_tot; i++) { - inputs.H_s.yx[k][j][i] = - PSTD.ca[n][i] * inputs.H_s.yx[k][j][i] + + inputs.H_s.yx(i, j, k) = + PSTD.ca[n][i] * inputs.H_s.yx(i, j, k) + PSTD.cb[n][i] * eh_vec[n][i][0] / ((double) PSTD.N_hx); } } @@ -2083,21 +2087,21 @@ void SimulationManager::execute() { k_loc = inputs.params.pml.Dzl + 1; } if (!inputs.materials[k][j][i]) { - inputs.H_s.yz[k][j][i] = - inputs.D.a.z[k_loc] * inputs.H_s.yz[k][j][i] + - inputs.D.b.z[k_loc] * (inputs.E_s.xy[k][j][i] + - inputs.E_s.xz[k][j][i] - - inputs.E_s.xy[k + 1][j][i] - - inputs.E_s.xz[k + 1][j][i]); + inputs.H_s.yz(i, j, k) = + inputs.D.a.z[k_loc] * inputs.H_s.yz(i, j, k) + + inputs.D.b.z[k_loc] * (inputs.E_s.xy(i, j, k) + + inputs.E_s.xz(i, j, k) - + inputs.E_s.xy(i, j, k + 1) - + inputs.E_s.xz(i, j, k + 1)); } else { - inputs.H_s.yz[k][j][i] = + inputs.H_s.yz(i, j, k) = inputs.Dmaterial.a.z[inputs.materials[k][j][i] - 1] * - inputs.H_s.yz[k][j][i] + + inputs.H_s.yz(i, j, k) + inputs.Dmaterial.b.z[inputs.materials[k][j][i] - 1] * - (inputs.E_s.xy[k][j][i] + - inputs.E_s.xz[k][j][i] - - inputs.E_s.xy[k + 1][j][i] - - inputs.E_s.xz[k + 1][j][i]); + (inputs.E_s.xy(i, j, k) + + inputs.E_s.xz(i, j, k) - + inputs.E_s.xy(i, j, k + 1) - + inputs.E_s.xz(i, j, k + 1)); } } } @@ -2130,33 +2134,33 @@ void SimulationManager::execute() { if (!inputs.materials[k][j][i]) { PSTD.ca[n][k] = inputs.D.a.z[k_loc]; PSTD.cb[n][k] = inputs.D.b.z[k_loc]; - // H_s.yz[k][j][i] = - // D.a.z[k_loc]*H_s.yz[k][j][i]+D.b.z[k_loc]*(E_s.xy[k][j][i] - // + E_s.xz[k][j][i] - E_s.xy[k+1][j][i] - E_s.xz[k+1][j][i]); + // H_s.yz(i,j,k) = + // D.a.z[k_loc]*H_s.yz(i,j,k)+D.b.z[k_loc]*(E_s.xy(i,j,k) + // + E_s.xz(i, j, k) - E_s.xy[k+1][j][i] - E_s.xz[k+1][j][i]); } else { PSTD.ca[n][k] = inputs.Dmaterial.a.z[inputs.materials[k][j][i] - 1]; PSTD.cb[n][k] = inputs.Dmaterial.b.z[inputs.materials[k][j][i] - 1]; - // H_s.yz[k][j][i] = - // Dmaterial.Da.z[materials[k][j][i]-1]*H_s.yz[k][j][i]+Dmaterial.Db.z[materials[k][j][i]-1]*(E_s.xy[k][j][i] - // + E_s.xz[k][j][i] - E_s.xy[k+1][j][i] - E_s.xz[k+1][j][i]); + // H_s.yz(i,j,k) = + // Dmaterial.Da.z[materials[k][j][i]-1]*H_s.yz(i,j,k)+Dmaterial.Db.z[materials[k][j][i]-1]*(E_s.xy(i,j,k) + // + E_s.xz(i, j, k) - E_s.xy[k+1][j][i] - E_s.xz[k+1][j][i]); } eh_vec[n][k][0] = - inputs.E_s.xy[k][j][i] + inputs.E_s.xz[k][j][i]; + inputs.E_s.xy(i, j, k) + inputs.E_s.xz(i, j, k); eh_vec[n][k][1] = 0.; } k = K_tot; - eh_vec[n][k][0] = inputs.E_s.xy[k][j][i] + inputs.E_s.xz[k][j][i]; + eh_vec[n][k][0] = inputs.E_s.xy(i, j, k) + inputs.E_s.xz(i, j, k); eh_vec[n][k][1] = 0.; first_derivative(eh_vec[n], eh_vec[n], PSTD.dk_hz, PSTD.N_hz, inputs.H_s.yz.plan_f[n], inputs.H_s.yz.plan_b[n]); for (k = 0; k < K_tot; k++) { - inputs.H_s.yz[k][j][i] = - PSTD.ca[n][k] * inputs.H_s.yz[k][j][i] - + inputs.H_s.yz(i, j, k) = + PSTD.ca[n][k] * inputs.H_s.yz(i, j, k) - PSTD.cb[n][k] * eh_vec[n][k][0] / ((double) PSTD.N_hz); } } @@ -2170,9 +2174,9 @@ void SimulationManager::execute() { for (k = 0; k <= K_tot; k++) for (j = 0; j < J_tot; j++) for (i = 0; i < (I_tot + 1); i++) - if (!inputs.materials[k][j][i]) inputs.H_s.xz[k][j][i] = 0.; + if (!inputs.materials[k][j][i]) inputs.H_s.xz(i, j, k) = 0.; else - inputs.H_s.xz[k][j][i] = 0.; + inputs.H_s.xz(i, j, k) = 0.; #pragma omp for // H_s.xy update @@ -2201,21 +2205,21 @@ void SimulationManager::execute() { else array_ind = (J_tot + 1) * k_loc + j; if (!inputs.materials[k][j][i]) - inputs.H_s.xy[k][j][i] = - inputs.D.a.y[array_ind] * inputs.H_s.xy[k][j][i] + - inputs.D.b.y[array_ind] * (inputs.E_s.zy[k][j][i] + - inputs.E_s.zx[k][j][i] - - inputs.E_s.zy[k][j + 1][i] - - inputs.E_s.zx[k][j + 1][i]); + inputs.H_s.xy(i, j, k) = + inputs.D.a.y[array_ind] * inputs.H_s.xy(i, j, k) + + inputs.D.b.y[array_ind] * (inputs.E_s.zy(i, j, k) + + inputs.E_s.zx(i, j, k) - + inputs.E_s.zy(i, j + 1, k) - + inputs.E_s.zx(i, j + 1, k)); else - inputs.H_s.xy[k][j][i] = + inputs.H_s.xy(i, j, k) = inputs.Dmaterial.a.y[inputs.materials[k][j][i] - 1] * - inputs.H_s.xy[k][j][i] + + inputs.H_s.xy(i, j, k) + inputs.Dmaterial.b.y[inputs.materials[k][j][i] - 1] * - (inputs.E_s.zy[k][j][i] + - inputs.E_s.zx[k][j][i] - - inputs.E_s.zy[k][j + 1][i] - - inputs.E_s.zx[k][j + 1][i]); + (inputs.E_s.zy(i, j, k) + + inputs.E_s.zx(i, j, k) - + inputs.E_s.zy(i, j + 1, k) - + inputs.E_s.zx(i, j + 1, k)); } #pragma omp for @@ -2245,30 +2249,30 @@ void SimulationManager::execute() { else array_ind = (I_tot + 1) * k_loc + i; if (!inputs.materials[k][j][i]) - inputs.H_s.yx[k][j][i] = - inputs.D.a.x[array_ind] * inputs.H_s.yx[k][j][i] + - inputs.D.b.x[array_ind] * (inputs.E_s.zx[k][j][i + 1] + - inputs.E_s.zy[k][j][i + 1] - - inputs.E_s.zx[k][j][i] - - inputs.E_s.zy[k][j][i]); + inputs.H_s.yx(i, j, k) = + inputs.D.a.x[array_ind] * inputs.H_s.yx(i, j, k) + + inputs.D.b.x[array_ind] * (inputs.E_s.zx(i + 1, j, k) + + inputs.E_s.zy(i + 1, j, k) - + inputs.E_s.zx(i, j, k) - + inputs.E_s.zy(i, j, k)); else - inputs.H_s.yx[k][j][i] = + inputs.H_s.yx(i, j, k) = inputs.Dmaterial.a.x[inputs.materials[k][j][i] - 1] * - inputs.H_s.yx[k][j][i] + + inputs.H_s.yx(i, j, k) + inputs.Dmaterial.b.x[inputs.materials[k][j][i] - 1] * - (inputs.E_s.zx[k][j][i + 1] + - inputs.E_s.zy[k][j][i + 1] - - inputs.E_s.zx[k][j][i] - - inputs.E_s.zy[k][j][i]); + (inputs.E_s.zx(i + 1, j, k) + + inputs.E_s.zy(i + 1, j, k) - + inputs.E_s.zx(i, j, k) - + inputs.E_s.zy(i, j, k)); } #pragma omp for for (k = 0; k <= K_tot; k++) { for (j = 0; j < (J_tot + 1); j++) for (i = 0; i < I_tot; i++) - if (!inputs.materials[k][j][i]) inputs.H_s.yz[k][j][i] = 0.; + if (!inputs.materials[k][j][i]) inputs.H_s.yz(i, j, k) = 0.; else - inputs.H_s.yz[k][j][i] = 0.; + inputs.H_s.yz(i, j, k) = 0.; } } @@ -2303,22 +2307,22 @@ void SimulationManager::execute() { else array_ind = (J_tot + 1) * k_loc + j; if (!inputs.materials[k][j][i]) - inputs.H_s.zy[k][j][i] = - inputs.D.a.y[array_ind] * inputs.H_s.zy[k][j][i] + + inputs.H_s.zy(i, j, k) = + inputs.D.a.y[array_ind] * inputs.H_s.zy(i, j, k) + inputs.D.b.y[array_ind] * - (inputs.E_s.xy[k][j + 1][i] + - inputs.E_s.xz[k][j + 1][i] - - inputs.E_s.xy[k][j][i] - - inputs.E_s.xz[k][j][i]); + (inputs.E_s.xy(i, j + 1, k) + + inputs.E_s.xz(i, j + 1, k) - + inputs.E_s.xy(i, j, k) - + inputs.E_s.xz(i, j, k)); else - inputs.H_s.zy[k][j][i] = + inputs.H_s.zy(i, j, k) = inputs.Dmaterial.a.y[inputs.materials[k][j][i] - 1] * - inputs.H_s.zy[k][j][i] + + inputs.H_s.zy(i, j, k) + inputs.Dmaterial.b.y[inputs.materials[k][j][i] - 1] * - (inputs.E_s.xy[k][j + 1][i] + - inputs.E_s.xz[k][j + 1][i] - - inputs.E_s.xy[k][j][i] - - inputs.E_s.xz[k][j][i]); + (inputs.E_s.xy(i, j + 1, k) + + inputs.E_s.xz(i, j + 1, k) - + inputs.E_s.xy(i, j, k) - + inputs.E_s.xz(i, j, k)); } // FDTD, H_s.zy } else { @@ -2351,25 +2355,25 @@ void SimulationManager::execute() { if (!inputs.materials[k][j][i]) { PSTD.ca[n][j] = inputs.D.a.y[array_ind]; PSTD.cb[n][j] = inputs.D.b.y[array_ind]; - // H_s.zy[k][j][i] = - // D.a.y[array_ind]*H_s.zy[k][j][i]+D.b.y[array_ind]*(E_s.xy[k][j+1][i] - // + E_s.xz[k][j+1][i] - E_s.xy[k][j][i] - E_s.xz[k][j][i]); + // H_s.zy(i,j,k) = + // D.a.y[array_ind]*H_s.zy(i,j,k)+D.b.y[array_ind]*(E_s.xy[k][j+1][i] + // + E_s.xz[k][j+1][i] - E_s.xy(i,j,k) - E_s.xz(i, j, k)); } else { PSTD.ca[n][j] = inputs.Dmaterial.a.y[inputs.materials[k][j][i] - 1]; PSTD.cb[n][j] = inputs.Dmaterial.b.y[inputs.materials[k][j][i] - 1]; - // H_s.zy[k][j][i] = - // Dmaterial.Da.y[materials[k][j][i]-1]*H_s.zy[k][j][i]+Dmaterial.Db.y[materials[k][j][i]-1]*(E_s.xy[k][j+1][i] - //+ E_s.xz[k][j+1][i] - E_s.xy[k][j][i] - E_s.xz[k][j][i]); + // H_s.zy(i,j,k) = + // Dmaterial.Da.y[materials[k][j][i]-1]*H_s.zy(i,j,k)+Dmaterial.Db.y[materials[k][j][i]-1]*(E_s.xy[k][j+1][i] + //+ E_s.xz[k][j+1][i] - E_s.xy(i,j,k) - E_s.xz(i, j, k)); } eh_vec[n][j][0] = - inputs.E_s.xy[k][j][i] + inputs.E_s.xz[k][j][i]; + inputs.E_s.xy(i, j, k) + inputs.E_s.xz(i, j, k); eh_vec[n][j][1] = 0.; } j = J_tot; - eh_vec[n][j][0] = inputs.E_s.xy[k][j][i] + inputs.E_s.xz[k][j][i]; + eh_vec[n][j][0] = inputs.E_s.xy(i, j, k) + inputs.E_s.xz(i, j, k); eh_vec[n][j][1] = 0.; first_derivative(eh_vec[n], eh_vec[n], PSTD.dk_hy, PSTD.N_hy, @@ -2377,8 +2381,8 @@ void SimulationManager::execute() { inputs.H_s.zy.plan_b[n]); for (j = 0; j < J_tot; j++) { - inputs.H_s.zy[k][j][i] = - PSTD.ca[n][j] * inputs.H_s.zy[k][j][i] + + inputs.H_s.zy(i, j, k) = + PSTD.ca[n][j] * inputs.H_s.zy(i, j, k) + PSTD.cb[n][j] * eh_vec[n][j][0] / ((double) PSTD.N_hy); } } @@ -2416,22 +2420,22 @@ void SimulationManager::execute() { else array_ind = (I_tot + 1) * k_loc + i; if (!inputs.materials[k][j][i]) - inputs.H_s.zx[k][j][i] = - inputs.D.a.x[array_ind] * inputs.H_s.zx[k][j][i] + + inputs.H_s.zx(i, j, k) = + inputs.D.a.x[array_ind] * inputs.H_s.zx(i, j, k) + inputs.D.b.x[array_ind] * - (inputs.E_s.yx[k][j][i] + - inputs.E_s.yz[k][j][i] - - inputs.E_s.yx[k][j][i + 1] - - inputs.E_s.yz[k][j][i + 1]); + (inputs.E_s.yx(i, j, k) + + inputs.E_s.yz(i, j, k) - + inputs.E_s.yx(i + 1, j, k) - + inputs.E_s.yz(i + 1, j, k)); else - inputs.H_s.zx[k][j][i] = + inputs.H_s.zx(i, j, k) = inputs.Dmaterial.a.x[inputs.materials[k][j][i] - 1] * - inputs.H_s.zx[k][j][i] + + inputs.H_s.zx(i, j, k) + inputs.Dmaterial.b.x[inputs.materials[k][j][i] - 1] * - (inputs.E_s.yx[k][j][i] + - inputs.E_s.yz[k][j][i] - - inputs.E_s.yx[k][j][i + 1] - - inputs.E_s.yz[k][j][i + 1]); + (inputs.E_s.yx(i, j, k) + + inputs.E_s.yz(i, j, k) - + inputs.E_s.yx(i + 1, j, k) - + inputs.E_s.yz(i + 1, j, k)); } // FDTD, H_s.zx } else { @@ -2462,15 +2466,15 @@ void SimulationManager::execute() { else array_ind = (I_tot + 1) * k_loc + i; if (!inputs.materials[k][j][i]) { - // H_s.zx[k][j][i] = - // D.a.x[array_ind]*H_s.zx[k][j][i]+D.b.x[array_ind]*(E_s.yx[k][j][i] - // + E_s.yz[k][j][i] - E_s.yx[k][j][i+1] - E_s.yz[k][j][i+1]); + // H_s.zx(i,j,k) = + // D.a.x[array_ind]*H_s.zx(i,j,k)+D.b.x[array_ind]*(E_s.yx(i,j,k) + // + E_s.yz(i,j,k) - E_s.yx[k][j][i+1] - E_s.yz[k][j][i+1]); PSTD.ca[n][i] = inputs.D.a.x[array_ind]; PSTD.cb[n][i] = inputs.D.b.x[array_ind]; } else { - // H_s.zx[k][j][i] = - // Dmaterial.Da.x[materials[k][j][i]-1]*H_s.zx[k][j][i]+Dmaterial.Db.x[materials[k][j][i]-1]*(E_s.yx[k][j][i] - //+ E_s.yz[k][j][i] - E_s.yx[k][j][i+1] - E_s.yz[k][j][i+1]); + // H_s.zx(i,j,k) = + // Dmaterial.Da.x[materials[k][j][i]-1]*H_s.zx(i,j,k)+Dmaterial.Db.x[materials[k][j][i]-1]*(E_s.yx(i,j,k) + //+ E_s.yz(i,j,k) - E_s.yx[k][j][i+1] - E_s.yz[k][j][i+1]); PSTD.ca[n][i] = inputs.Dmaterial.a.x[inputs.materials[k][j][i] - 1]; PSTD.cb[n][i] = @@ -2478,11 +2482,11 @@ void SimulationManager::execute() { } eh_vec[n][i][0] = - inputs.E_s.yx[k][j][i] + inputs.E_s.yz[k][j][i]; + inputs.E_s.yx(i, j, k) + inputs.E_s.yz(i, j, k); eh_vec[n][i][1] = 0.; } i = I_tot; - eh_vec[n][i][0] = inputs.E_s.yx[k][j][i] + inputs.E_s.yz[k][j][i]; + eh_vec[n][i][0] = inputs.E_s.yx(i, j, k) + inputs.E_s.yz(i, j, k); eh_vec[n][i][1] = 0.; @@ -2491,8 +2495,8 @@ void SimulationManager::execute() { inputs.H_s.zx.plan_b[n]); for (i = 0; i < I_tot; i++) { - inputs.H_s.zx[k][j][i] = - PSTD.ca[n][i] * inputs.H_s.zx[k][j][i] - + inputs.H_s.zx(i, j, k) = + PSTD.ca[n][i] * inputs.H_s.zx(i, j, k) - PSTD.cb[n][i] * eh_vec[n][i][0] / ((double) PSTD.N_hx); } } diff --git a/tdms/src/simulation_manager/execute_update_Ex_split.cpp b/tdms/src/simulation_manager/execute_update_Ex_split.cpp index eec16cd94..57d6f6e26 100644 --- a/tdms/src/simulation_manager/execute_update_Ex_split.cpp +++ b/tdms/src/simulation_manager/execute_update_Ex_split.cpp @@ -120,61 +120,61 @@ void SimulationManager::update_Exy(LoopVariables &lv) { double Enp1, Jnp1; if (solver_method == SolverMethod::FiniteDifference) { - Enp1 = Ca * inputs.E_s.xy[k][j][i] + - Cb * (inputs.H_s.zy[k][j][i] + inputs.H_s.zx[k][j][i] - - inputs.H_s.zy[k][j - 1][i] - inputs.H_s.zx[k][j - 1][i]); + Enp1 = Ca * inputs.E_s.xy(i, j, k) + + Cb * (inputs.H_s.zy(i, j, k) + inputs.H_s.zx(i, j, k) - + inputs.H_s.zy(i, j - 1, k) - inputs.H_s.zx(i, j - 1, k)); if ((lv.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * lv.E_nm1.xy[k][j][i] - + Enp1 += Cc * lv.E_nm1.xy(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dy * - ((1 + alpha_l) * lv.J_s.xy[k][j][i] + - beta_l * lv.J_nm1.xy[k][j][i]); + ((1 + alpha_l) * lv.J_s.xy(i, j, k) + + beta_l * lv.J_nm1.xy(i, j, k)); if (lv.is_conductive && rho) - Enp1 += Cb * inputs.params.delta.dy * lv.J_c.xy[k][j][i]; + Enp1 += Cb * inputs.params.delta.dy * lv.J_c.xy(i, j, k); if ((lv.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * lv.J_s.xy[k][j][i] + - beta_l * lv.J_nm1.xy[k][j][i] + + Jnp1 = alpha_l * lv.J_s.xy(i, j, k) + + beta_l * lv.J_nm1.xy(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - lv.E_nm1.xy[k][j][i]); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.xy[k][j][i]; + (Enp1 - lv.E_nm1.xy(i, j, k)); + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.xy(i, j, k); - lv.E_nm1.xy[k][j][i] = inputs.E_s.xy[k][j][i]; - lv.J_nm1.xy[k][j][i] = lv.J_s.xy[k][j][i]; - lv.J_s.xy[k][j][i] = Jnp1; + lv.E_nm1.xy(i, j, k) = inputs.E_s.xy(i, j, k); + lv.J_nm1.xy(i, j, k) = lv.J_s.xy(i, j, k); + lv.J_s.xy(i, j, k) = Jnp1; } if (lv.is_conductive && rho) { - lv.J_c.xy[k][j][i] -= rho * (Enp1 + inputs.E_s.xy[k][j][i]); + lv.J_c.xy(i, j, k) -= rho * (Enp1 + inputs.E_s.xy(i, j, k)); } - inputs.E_s.xy[k][j][i] = Enp1; + inputs.E_s.xy(i, j, k) = Enp1; } else {// pseudo-spectral Enp1 = 0.0; - // Enp1 = Ca*E_s.xy[k][j][i]+Cb*(H_s.zy[k][j][i] + H_s.zx[k][j][i] - + // Enp1 = Ca*E_s.xy(i,j,k)+Cb*(H_s.zy(i,j,k) + H_s.zx(i,j,k) - // H_s.zy[k][j-1][i] - H_s.zx[k][j-1][i]); if ((lv.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * lv.E_nm1.xy[k][j][i] - + Enp1 += Cc * lv.E_nm1.xy(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dy * - ((1 + alpha_l) * lv.J_s.xy[k][j][i] + - beta_l * lv.J_nm1.xy[k][j][i]); + ((1 + alpha_l) * lv.J_s.xy(i, j, k) + + beta_l * lv.J_nm1.xy(i, j, k)); if (lv.is_conductive && rho) - Enp1 += Cb * inputs.params.delta.dy * lv.J_c.xy[k][j][i]; + Enp1 += Cb * inputs.params.delta.dy * lv.J_c.xy(i, j, k); if ((lv.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * lv.J_s.xy[k][j][i] + - beta_l * lv.J_nm1.xy[k][j][i] + + Jnp1 = alpha_l * lv.J_s.xy(i, j, k) + + beta_l * lv.J_nm1.xy(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - lv.E_nm1.xy[k][j][i]); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.xy[k][j][i]; + (Enp1 - lv.E_nm1.xy(i, j, k)); + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.xy(i, j, k); - lv.E_nm1.xy[k][j][i] = inputs.E_s.xy[k][j][i]; - lv.J_nm1.xy[k][j][i] = lv.J_s.xy[k][j][i]; - lv.J_s.xy[k][j][i] = Jnp1; + lv.E_nm1.xy(i, j, k) = inputs.E_s.xy(i, j, k); + lv.J_nm1.xy(i, j, k) = lv.J_s.xy(i, j, k); + lv.J_s.xy(i, j, k) = Jnp1; } if (lv.is_conductive && rho) { - lv.J_c.xy[k][j][i] -= rho * (Enp1 + inputs.E_s.xy[k][j][i]); + lv.J_c.xy(i, j, k) -= rho * (Enp1 + inputs.E_s.xy(i, j, k)); } - eh_vec[n][j][0] = inputs.H_s.zy[k][j][i] + inputs.H_s.zx[k][j][i]; + eh_vec[n][j][0] = inputs.H_s.zy(i, j, k) + inputs.H_s.zx(i, j, k); eh_vec[n][j][1] = 0.; PSTD.ca[n][j - 1] = Ca; PSTD.cb[n][j - 1] = Cb; @@ -182,13 +182,13 @@ void SimulationManager::update_Exy(LoopVariables &lv) { } if (solver_method == SolverMethod::PseudoSpectral && J_tot > 1) { int j = 0; - eh_vec[n][j][0] = inputs.H_s.zy[k][j][i] + inputs.H_s.zx[k][j][i]; + eh_vec[n][j][0] = inputs.H_s.zy(i, j, k) + inputs.H_s.zx(i, j, k); eh_vec[n][j][1] = 0.; first_derivative(eh_vec[n], eh_vec[n], PSTD.dk_ey, PSTD.N_ey, inputs.E_s.xy.plan_f[n], inputs.E_s.xy.plan_b[n]); for (j = 1; j < J_tot; j++) { - inputs.E_s.xy[k][j][i] = - PSTD.ca[n][j - 1] * inputs.E_s.xy[k][j][i] + + inputs.E_s.xy(i, j, k) = + PSTD.ca[n][j - 1] * inputs.E_s.xy(i, j, k) + PSTD.cb[n][j - 1] * eh_vec[n][j][0] / ((double) PSTD.N_ey); } } @@ -302,58 +302,58 @@ void SimulationManager::update_Exz(LoopVariables &lv) { // variables... double Enp1, Jnp1; if (solver_method == SolverMethod::FiniteDifference) { - Enp1 = Ca * inputs.E_s.xz[k][j][i] + - Cb * (inputs.H_s.yx[k - 1][j][i] + inputs.H_s.yz[k - 1][j][i] - - inputs.H_s.yx[k][j][i] - inputs.H_s.yz[k][j][i]); + Enp1 = Ca * inputs.E_s.xz(i, j, k) + + Cb * (inputs.H_s.yx(i, j, k - 1) + inputs.H_s.yz(i, j, k - 1) - + inputs.H_s.yx(i, j, k) - inputs.H_s.yz(i, j, k)); if ((lv.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * lv.E_nm1.xz[k][j][i] - + Enp1 += Cc * lv.E_nm1.xz(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dz * - ((1 + alpha_l) * lv.J_s.xz[k][j][i] + - beta_l * lv.J_nm1.xz[k][j][i]); + ((1 + alpha_l) * lv.J_s.xz(i, j, k) + + beta_l * lv.J_nm1.xz(i, j, k)); if (lv.is_conductive && rho) - Enp1 += Cb * inputs.params.delta.dz * lv.J_c.xz[k][j][i]; + Enp1 += Cb * inputs.params.delta.dz * lv.J_c.xz(i, j, k); if ((lv.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * lv.J_s.xz[k][j][i] + - beta_l * lv.J_nm1.xz[k][j][i] + + Jnp1 = alpha_l * lv.J_s.xz(i, j, k) + + beta_l * lv.J_nm1.xz(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - lv.E_nm1.xz[k][j][i]); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.xz[k][j][i]; - lv.E_nm1.xz[k][j][i] = inputs.E_s.xz[k][j][i]; - lv.J_nm1.xz[k][j][i] = lv.J_s.xz[k][j][i]; - lv.J_s.xz[k][j][i] = Jnp1; + (Enp1 - lv.E_nm1.xz(i, j, k)); + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.xz(i, j, k); + lv.E_nm1.xz(i, j, k) = inputs.E_s.xz(i, j, k); + lv.J_nm1.xz(i, j, k) = lv.J_s.xz(i, j, k); + lv.J_s.xz(i, j, k) = Jnp1; } if (lv.is_conductive && rho) { - lv.J_c.xz[k][j][i] -= rho * (Enp1 + inputs.E_s.xz[k][j][i]); + lv.J_c.xz(i, j, k) -= rho * (Enp1 + inputs.E_s.xz(i, j, k)); } - inputs.E_s.xz[k][j][i] = Enp1; + inputs.E_s.xz(i, j, k) = Enp1; } else {// psuedo-spectral - // Enp1 = Ca*E_s.xz[k][j][i]+Cb*(H_s.yx[k-1][j][i] + H_s.yz[k-1][j][i] - // - H_s.yx[k][j][i] - H_s.yz[k][j][i]); + // Enp1 = Ca*E_s.xz(i,j,k)+Cb*(H_s.yx[k-1][j][i] + H_s.yz[k-1][j][i] + // - H_s.yx(i,j,k) - H_s.yz(i,j,k)); if ((lv.is_dispersive || inputs.params.is_disp_ml) && gamma_l) - Enp1 += Cc * lv.E_nm1.xz[k][j][i] - + Enp1 += Cc * lv.E_nm1.xz(i, j, k) - 1. / 2. * Cb * inputs.params.delta.dz * - ((1 + alpha_l) * lv.J_s.xz[k][j][i] + - beta_l * lv.J_nm1.xz[k][j][i]); + ((1 + alpha_l) * lv.J_s.xz(i, j, k) + + beta_l * lv.J_nm1.xz(i, j, k)); if (lv.is_conductive && rho) - Enp1 += Cb * inputs.params.delta.dz * lv.J_c.xz[k][j][i]; + Enp1 += Cb * inputs.params.delta.dz * lv.J_c.xz(i, j, k); if ((lv.is_dispersive || inputs.params.is_disp_ml) && gamma_l) { - Jnp1 = alpha_l * lv.J_s.xz[k][j][i] + - beta_l * lv.J_nm1.xz[k][j][i] + + Jnp1 = alpha_l * lv.J_s.xz(i, j, k) + + beta_l * lv.J_nm1.xz(i, j, k) + kappa_l * gamma_l / (2. * inputs.params.dt) * - (Enp1 - lv.E_nm1.xz[k][j][i]); - Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.xz[k][j][i]; - lv.E_nm1.xz[k][j][i] = inputs.E_s.xz[k][j][i]; - lv.J_nm1.xz[k][j][i] = lv.J_s.xz[k][j][i]; - lv.J_s.xz[k][j][i] = Jnp1; + (Enp1 - lv.E_nm1.xz(i, j, k)); + Jnp1 += sigma_l / EPSILON0 * gamma_l * inputs.E_s.xz(i, j, k); + lv.E_nm1.xz(i, j, k) = inputs.E_s.xz(i, j, k); + lv.J_nm1.xz(i, j, k) = lv.J_s.xz(i, j, k); + lv.J_s.xz(i, j, k) = Jnp1; } if (lv.is_conductive && rho) { - lv.J_c.xz[k][j][i] -= rho * (Enp1 + inputs.E_s.xz[k][j][i]); + lv.J_c.xz(i, j, k) -= rho * (Enp1 + inputs.E_s.xz(i, j, k)); } - eh_vec[n][k][0] = inputs.H_s.yx[k][j][i] + inputs.H_s.yz[k][j][i]; + eh_vec[n][k][0] = inputs.H_s.yx(i, j, k) + inputs.H_s.yz(i, j, k); eh_vec[n][k][1] = 0.; PSTD.ca[n][k - 1] = Ca; PSTD.cb[n][k - 1] = Cb; @@ -361,15 +361,15 @@ void SimulationManager::update_Exz(LoopVariables &lv) { } if (solver_method == SolverMethod::PseudoSpectral) { int k = 0; - eh_vec[n][k][0] = inputs.H_s.yx[k][j][i] + inputs.H_s.yz[k][j][i]; + eh_vec[n][k][0] = inputs.H_s.yx(i, j, k) + inputs.H_s.yz(i, j, k); eh_vec[n][k][1] = 0.; first_derivative(eh_vec[n], eh_vec[n], PSTD.dk_ez, PSTD.N_ez, inputs.E_s.xz.plan_f[n], inputs.E_s.xz.plan_b[n]); for (k = 1; k < K_tot; k++) { - inputs.E_s.xz[k][j][i] = - PSTD.ca[n][k - 1] * inputs.E_s.xz[k][j][i] - + inputs.E_s.xz(i, j, k) = + PSTD.ca[n][k - 1] * inputs.E_s.xz(i, j, k) - PSTD.cb[n][k - 1] * eh_vec[n][k][0] / ((double) PSTD.N_ez); } } diff --git a/tdms/src/simulation_manager/source_updates_pulsed.cpp b/tdms/src/simulation_manager/source_updates_pulsed.cpp index af1742b06..dcbe41136 100644 --- a/tdms/src/simulation_manager/source_updates_pulsed.cpp +++ b/tdms/src/simulation_manager/source_updates_pulsed.cpp @@ -139,7 +139,7 @@ void SimulationManager::update_source_terms_pulsed(double time_E, int tind) { d_constant * common_amplitude * real(source_value * common_phase); // Update broadband source term if (inputs.params.eyi_present) { - inputs.H_s.xz[cell_to_update] -= d_constant * inputs.Ei.y[tind][j][i]; + inputs.H_s.xz[cell_to_update] -= d_constant * inputs.Ei.y(i, j, tind); } } for (int i = 0; i < I_tot; i++) { @@ -152,7 +152,7 @@ void SimulationManager::update_source_terms_pulsed(double time_E, int tind) { d_constant * common_amplitude * real(source_value * common_phase); // Update broadband source term if (inputs.params.exi_present) { - inputs.H_s.yz[cell_to_update] += d_constant * inputs.Ei.x[tind][j][i]; + inputs.H_s.yz[cell_to_update] += d_constant * inputs.Ei.x(i, j, tind); } } } else { @@ -165,7 +165,7 @@ void SimulationManager::update_source_terms_pulsed(double time_E, int tind) { inputs.H_s.xz[cell_to_update] -= d_constant * common_amplitude * real(source_value * common_phase); if (inputs.params.eyi_present) { - inputs.H_s.xz[cell_to_update] -= d_constant * inputs.Ei.y[tind][j][i]; + inputs.H_s.xz[cell_to_update] -= d_constant * inputs.Ei.y(i, j, tind); } } } @@ -178,7 +178,7 @@ void SimulationManager::update_source_terms_pulsed(double time_E, int tind) { inputs.H_s.yz[cell_to_update] += d_constant * common_amplitude * real(source_value * common_phase); if (inputs.params.exi_present) { - inputs.H_s.yz[cell_to_update] += d_constant * inputs.Ei.x[tind][j][i]; + inputs.H_s.yz[cell_to_update] += d_constant * inputs.Ei.x(i, j, tind); } } } diff --git a/tdms/src/simulation_parameters.cpp b/tdms/src/simulation_parameters.cpp index 26a0857a4..33803f83c 100644 --- a/tdms/src/simulation_parameters.cpp +++ b/tdms/src/simulation_parameters.cpp @@ -5,6 +5,8 @@ #include +#include "arrays/incident_field.h" + using namespace std; SimulationParameters::SimulationParameters() = default; @@ -126,7 +128,8 @@ void SimulationParameters::unpack_from_input_matrices( exi_present = Ei.x.has_elements(); eyi_present = Ei.y.has_elements(); - // set_Np + // set_Np -> do we actually USE this? f_ex_vec goes out of scope immediately + // after... FrequencyExtractVector f_ex_vec(in_matrices["f_ex_vec"], omega_an); set_Np(f_ex_vec); } diff --git a/tdms/tests/unit/array_tests/test_DTilde.cpp b/tdms/tests/unit/array_tests/test_DTilde.cpp index edf804fd7..e38f4d589 100644 --- a/tdms/tests/unit/array_tests/test_DTilde.cpp +++ b/tdms/tests/unit/array_tests/test_DTilde.cpp @@ -8,6 +8,7 @@ #include "array_test_class.h" #include "arrays.h" +#include "arrays/dtilde.h" using namespace std; diff --git a/tdms/tests/unit/array_tests/test_IncidentField.cpp b/tdms/tests/unit/array_tests/test_IncidentField.cpp index c81722b36..03557e796 100644 --- a/tdms/tests/unit/array_tests/test_IncidentField.cpp +++ b/tdms/tests/unit/array_tests/test_IncidentField.cpp @@ -7,7 +7,7 @@ #include #include "array_test_class.h" -#include "arrays.h" +#include "arrays/incident_field.h" #include "unit_test_utils.h" using namespace std; @@ -73,9 +73,9 @@ void IncidentFieldTest::test_correct_construction() { for (int kk = 0; kk < n_layers; kk++) { elements_set_correctly = elements_set_correctly && - (abs(i_field.x[kk][jj][ii] - + (abs(i_field.x(ii, jj, kk) - 1. / ((double) (kk + jj + ii + 1))) < TOLERANCE) && - (abs(i_field.y[kk][jj][ii] - + (abs(i_field.y(ii, jj, kk) - 1. / ((double) (kk + jj + ii + 1))) < TOLERANCE); } } diff --git a/tdms/tests/unit/array_tests/test_Tensor3D.cpp b/tdms/tests/unit/array_tests/test_Tensor3D.cpp index a661b4d59..514691aed 100644 --- a/tdms/tests/unit/array_tests/test_Tensor3D.cpp +++ b/tdms/tests/unit/array_tests/test_Tensor3D.cpp @@ -15,39 +15,71 @@ using tdms_tests::TOLERANCE; void Tensor3DTest::test_correct_construction() { - Tensor3D *t3d; + Tensor3D t3d; SECTION("Default constructor") { // default constructor should assign all dimensions to 0, and the tensor - // itself should be a nullptr - t3d = new Tensor3D; - // should have no elements - REQUIRE(!(t3d->has_elements())); + // should have no elements as a result + REQUIRE(!(t3d.has_elements())); // let us assign some free memory to this tensor via allocate() - t3d->allocate(n_layers, n_cols, n_rows); + t3d.allocate(n_layers, n_cols, n_rows); // we should now "have elements", even though they are unassigned - REQUIRE(t3d->has_elements()); + REQUIRE(t3d.has_elements()); } SECTION("Overloaded constructor") { - // also implicitly tests initialise() - double ***p = (double ***) malloc(n_layers * sizeof(double **)); - for (int k = 0; k < n_layers; k++) { - p[k] = (double **) malloc(n_cols * sizeof(double *)); - for (int j = 0; j < n_cols; j++) { - p[k][j] = (double *) malloc(n_rows * sizeof(double)); + // testing initialise() as the overloaded constructor is a wrapper for this + int ***p = nullptr; + // assemble p such that the [i][j][k]th entry is the strided vector index, + // i*n_layers*n_cols + j*n_layers + k + SECTION("Reverse buffer") { + // Test when we read in a 3D buffer of shape [n_layers][n_cols][n_rows] + p = (int ***) malloc(n_layers * sizeof(int **)); + for (int k = 0; k < n_layers; k++) { + p[k] = (int **) malloc(n_cols * sizeof(int *)); + for (int j = 0; j < n_cols; j++) { + p[k][j] = (int *) malloc(n_rows * sizeof(int)); + for (int i = 0; i < n_rows; i++) { + p[k][j][i] = i * n_layers * n_cols + j * n_layers + k; + } + } + } + t3d.initialise(p, n_layers, n_cols, n_rows, true); + // Malloc'd values are copied, so we can tear down now + free(p); + } + SECTION("Ordered buffer") { + // Test when we read in a 3D buffer of shape [n_rows][n_cols][n_layers] + p = (int ***) malloc(n_rows * sizeof(int **)); + for (int i = 0; i < n_rows; i++) { + p[i] = (int **) malloc(n_cols * sizeof(int *)); + for (int j = 0; j < n_cols; j++) { + p[i][j] = (int *) malloc(n_rows * sizeof(int)); + for (int k = 0; k < n_layers; k++) { + p[i][j][k] = i * n_layers * n_cols + j * n_layers + k; + } + } } + t3d.initialise(p, n_layers, n_cols, n_rows, false); + // Malloc'd values are copied, so we can tear down now + free(p); } - t3d = new Tensor3D(p, n_layers, n_cols, n_rows); // this tensor should be flagged as "having elements", since we provided a // pointer in the constructor - REQUIRE(t3d->has_elements()); + bool elements_are_correct = t3d.has_elements(); + for (int i = 0; i < n_rows; i++) { + for (int j = 0; j < n_cols; j++) { + for (int k = 0; k < n_layers; k++) { + elements_are_correct = + elements_are_correct && + (t3d(i, j, k) == i * n_layers * n_cols + j * n_layers + k); + } + } + } + REQUIRE(t3d.has_elements()); } - // tear down assigned memory - delete t3d; } void Tensor3DTest::test_other_methods() { - Tensor3D t3d; - t3d.allocate(n_layers, n_cols, n_rows); + Tensor3D t3d(n_layers, n_cols, n_rows); t3d.zero(); SECTION("allocate() and zero()") { // we should be able to flag this tensor has elements, so the bool should be @@ -57,7 +89,7 @@ void Tensor3DTest::test_other_methods() { for (int j = 0; j < n_cols; j++) { for (int i = 0; i < n_rows; i++) { allocated_and_zero = - allocated_and_zero && (abs(t3d[k][j][i]) < TOLERANCE); + allocated_and_zero && (abs(t3d(i, j, k)) < TOLERANCE); } } } @@ -71,7 +103,7 @@ void Tensor3DTest::test_other_methods() { // 1, so the analytic norm is 16. for (int k = 0; k < n_layers; k++) { for (int j = 0; j < n_cols; j++) { - for (int i = 0; i < n_rows; i++) { t3d[k][j][i] = (i + j + k) % 2; } + for (int i = 0; i < n_rows; i++) { t3d(i, j, k) = (i + j + k) % 2; } } } // the analytic frobenuis norm of the tensor values diff --git a/tdms/tests/unit/field_tests/test_Field_interpolation.cpp b/tdms/tests/unit/field_tests/test_Field_interpolation.cpp index fe02de41c..94fb639c9 100644 --- a/tdms/tests/unit/field_tests/test_Field_interpolation.cpp +++ b/tdms/tests/unit/field_tests/test_Field_interpolation.cpp @@ -27,12 +27,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 @@ -78,7 +78,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 @@ -112,12 +112,12 @@ TEST_CASE("E-field interpolation check") { E.imag.z[kk][jj][ii] = 0.; // assign to "time domain" ElectricSplitField - use weighting that sums // to 1 to check addition is behaving as planned - E_split.xy[kk][jj][ii] = x_comp_value; - E_split.xz[kk][jj][ii] = 0.; - E_split.yx[kk][jj][ii] = y_comp_value * .5; - E_split.yz[kk][jj][ii] = y_comp_value * .5; - E_split.zx[kk][jj][ii] = z_comp_value * .25; - E_split.zy[kk][jj][ii] = z_comp_value * .75; + E_split.xy(ii, jj, kk) = x_comp_value; + E_split.xz(ii, jj, kk) = 0.; + E_split.yx(ii, jj, kk) = y_comp_value * .5; + E_split.yz(ii, jj, kk) = y_comp_value * .5; + E_split.zx(ii, jj, kk) = z_comp_value * .25; + E_split.zy(ii, jj, kk) = z_comp_value * .75; } } } @@ -131,14 +131,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) @@ -160,8 +158,8 @@ TEST_CASE("E-field interpolation check") { double Ex_interp = E.interpolate_to_centre_of(AxialDirection::X, current_cell) .real(); - Ex_error[kk][jj][ii - 1] = Ex_interp - Ex_exact; - Ex_split_error[kk][jj][ii - 1] = Ex_split_interp - Ex_exact; + Ex_error(ii - 1, jj, kk) = Ex_interp - Ex_exact; + Ex_split_error(ii - 1, jj, kk) = Ex_split_interp - Ex_exact; } // Ey interpolation @@ -172,8 +170,8 @@ TEST_CASE("E-field interpolation check") { double Ey_interp = E.interpolate_to_centre_of(AxialDirection::Y, current_cell) .real(); - Ey_error[kk][jj - 1][ii] = Ey_interp - Ey_exact; - Ey_split_error[kk][jj - 1][ii] = Ey_split_interp - Ey_exact; + Ey_error(ii, jj - 1, kk) = Ey_interp - Ey_exact; + Ey_split_error(ii, jj - 1, kk) = Ey_split_interp - Ey_exact; } // Ez interpolation @@ -184,8 +182,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; } } } @@ -207,8 +205,8 @@ TEST_CASE("E-field interpolation check") { // as such, make a new array for safety double jk_errors[Nx - 1], jk_split_errors[Nx - 1]; for (int ii = 0; ii < Nx - 1; ii++) { - jk_errors[ii] = Ex_error[kk][jj][ii]; - jk_split_errors[ii] = Ex_split_error[kk][jj][ii]; + jk_errors[ii] = Ex_error(ii, jj, kk); + jk_split_errors[ii] = Ex_split_error(ii, jj, kk); } // compute norm-error of this slice double jk_slice_error = euclidean(jk_errors, Nx - 1), @@ -225,8 +223,8 @@ TEST_CASE("E-field interpolation check") { for (int kk = 0; kk < Nz; kk++) { double ik_errors[Ny - 1], ik_split_errors[Ny - 1]; for (int jj = 0; jj < Ny - 1; jj++) { - ik_errors[jj] = Ey_error[kk][jj][ii]; - ik_split_errors[jj] = Ey_split_error[kk][jj][ii]; + ik_errors[jj] = Ey_error(ii, jj, kk); + ik_split_errors[jj] = Ey_split_error(ii, jj, kk); } double ik_slice_error = euclidean(ik_errors, Ny - 1), ik_split_slice_error = euclidean(ik_split_errors, Ny - 1); @@ -241,8 +239,8 @@ TEST_CASE("E-field interpolation check") { for (int jj = 0; jj < Ny; jj++) { double ij_errors[Nz - 1], ij_split_errors[Nz - 1]; for (int kk = 0; kk < Nz - 1; kk++) { - ij_errors[kk] = Ez_error[kk][jj][ii]; - ij_split_errors[kk] = Ez_split_error[kk][jj][ii]; + ij_errors[kk] = Ez_error(ii, jj, kk); + ij_split_errors[kk] = Ez_split_error(ii, jj, kk); } double ij_slice_error = euclidean(ij_errors, Nz - 1), ij_split_slice_error = euclidean(ij_split_errors, Nz - 1); @@ -295,7 +293,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 @@ -356,12 +354,12 @@ TEST_CASE("H-field interpolation check") { H.real.z[kk][jj][ii] = 0.; // assign to "time domain" ElectricSplitField - use weighting that sums // to 1 to check addition is behaving as planned - H_split.xy[kk][jj][ii] = x_comp_value; - H_split.xz[kk][jj][ii] = 0.; - H_split.yx[kk][jj][ii] = y_comp_value * .25; - H_split.yz[kk][jj][ii] = y_comp_value * .75; - H_split.zx[kk][jj][ii] = z_comp_value * .125; - H_split.zy[kk][jj][ii] = z_comp_value * .875; + H_split.xy(ii, jj, kk) = x_comp_value; + H_split.xz(ii, jj, kk) = 0.; + H_split.yx(ii, jj, kk) = y_comp_value * .25; + H_split.yz(ii, jj, kk) = y_comp_value * .75; + H_split.zx(ii, jj, kk) = z_comp_value * .125; + H_split.zy(ii, jj, kk) = z_comp_value * .875; } } } @@ -405,8 +403,8 @@ TEST_CASE("H-field interpolation check") { double Hx_interp = H.interpolate_to_centre_of(AxialDirection::X, current_cell) .imag(); - Hx_error[kk - 1][jj - 1][ii] = Hx_interp - Hx_exact; - Hx_split_error[kk - 1][jj - 1][ii] = Hx_split_interp - Hx_exact; + Hx_error(ii, jj - 1, kk - 1) = Hx_interp - Hx_exact; + Hx_split_error(ii, jj - 1, kk - 1) = Hx_split_interp - Hx_exact; } // Hy interpolation @@ -417,8 +415,8 @@ TEST_CASE("H-field interpolation check") { double Hy_interp = H.interpolate_to_centre_of(AxialDirection::Y, current_cell) .imag(); - Hy_error[kk - 1][jj][ii - 1] = Hy_interp - Hy_exact; - Hy_split_error[kk - 1][jj][ii - 1] = Hy_split_interp - Hy_exact; + Hy_error(ii - 1, jj, kk - 1) = Hy_interp - Hy_exact; + Hy_split_error(ii - 1, jj, kk - 1) = Hy_split_interp - Hy_exact; } // Hz interpolation @@ -429,8 +427,8 @@ TEST_CASE("H-field interpolation check") { double Hz_interp = H.interpolate_to_centre_of(AxialDirection::Z, current_cell) .imag(); - Hz_error[kk][jj - 1][ii - 1] = Hz_interp - Hz_exact; - Hz_split_error[kk][jj - 1][ii - 1] = Hz_split_interp - Hz_exact; + Hz_error(ii - 1, jj - 1, kk) = Hz_interp - Hz_exact; + Hz_split_error(ii - 1, jj - 1, kk) = Hz_split_interp - Hz_exact; } } } diff --git a/tdms/tests/unit/field_tests/test_SplitField.cpp b/tdms/tests/unit/field_tests/test_SplitField.cpp index fb069796b..9b6ebb454 100644 --- a/tdms/tests/unit/field_tests/test_SplitField.cpp +++ b/tdms/tests/unit/field_tests/test_SplitField.cpp @@ -29,9 +29,9 @@ TEST_CASE("SplitField: allocate and zero") { for (int k = 0; k < 1; k++) { for (int j = 0; j < 2; j++) { for (int i = 0; i < 3; i++) { - REQUIRE(is_close(field.xy[k][j][i], 0.0)); - REQUIRE(is_close(field.zy[k][j][i], 0.0)); - REQUIRE(is_close(field.yz[k][j][i], 0.0)); + REQUIRE(is_close(field.xy(i, j, k), 0.0)); + REQUIRE(is_close(field.zy(i, j, k), 0.0)); + REQUIRE(is_close(field.yz(i, j, k), 0.0)); } } } diff --git a/tdms/tests/unit/field_tests/test_SplitFieldComponent.cpp b/tdms/tests/unit/field_tests/test_SplitFieldComponent.cpp index 229c5e7e8..df7e3c1dd 100644 --- a/tdms/tests/unit/field_tests/test_SplitFieldComponent.cpp +++ b/tdms/tests/unit/field_tests/test_SplitFieldComponent.cpp @@ -19,8 +19,8 @@ TEST_CASE("SplitFieldComponent: set component") { tmp.allocate(2, 2, 2);// 2x2x2 tensor tmp.zero(); - REQUIRE(is_close(tmp[0][1][0], 0.0)); + REQUIRE(is_close(tmp(0, 1, 0), 0.0)); - tmp[0][1][0] = 1.0; - REQUIRE(is_close(tmp[0][1][0], 1.0)); + tmp(0, 1, 0) = 1.0; + REQUIRE(is_close(tmp(0, 1, 0), 1.0)); }