From a60f7be727e7dfb1e85838253d3447b34f40e145 Mon Sep 17 00:00:00 2001 From: Will Graham <32364977+willGraham01@users.noreply.github.com> Date: Fri, 16 Dec 2022 15:30:58 +0000 Subject: [PATCH] Replace cubic interpolation with Band-Limited Interpolation (#151) * Standardise the interpolation notation * Replace interpolation functions (#144) * Replace interpolation calls in omp loop L#1171-1193 * Successfully replace split-field interpolation (pending Mac brew fix) * Highlight fieldsample for extraction * Remove debug variables * Remove old interpolation call to bandlimited method * Use simpler notation in loop * replace 3D interpolation, iterator.cpp L4545 * Replace all calls on L4545 * Replace all interpolation calls * Reduce code bloat by moving up class hierachy * Remove old interpolation function files * Change extraction logic to avoid seg-faults * Reintroduce functionality, found error * Now everything seems to work * Interpolation is working! * Cleanup unused functions * Fix typo in MagneticSplitField's Hz-interpolation proceedure * Code formatting * Apply suggestions from code review Have avoided renaming functions, as that needs to be done consistently on my local branch. Co-authored-by: Sam Cunliffe * Update tdms/include/field.h Missed one Co-authored-by: Sam Cunliffe * Rename TE/TM functions to full-names * Random indent offset fix * TE/TM to full names (@samcunliffe) * Fix missing template * Adding a `CellCoordinate` class to group `i,j,k` (#149) * CellCoordinate into every declaration in field.h that has i,j,k * Field (and subclass) methods only take CellCoordinates now * Pragma once for cell_coordinate * Update tests to use cell_coordinate * CellCoordinate into every declaration in field.h that has i,j,k * Field (and subclass) methods only take CellCoordinates now * Pragma once for cell_coordinate * Update tests to use cell_coordinate * Catch stray j variable * Add docstrings to CellCoordinate * @samcunliffe's line shortening * Catch offset change in function definition * Remove additional logger statement * Interpolation method as a toggle (#199) * Add command line option to toggle interpolation method * Docment command-line option * Actually add ability to toggle behaviour * Update interpolation determination test * Sam's suggestions * Fix main test merge in Co-authored-by: willGraham01 <1willgraham@gmail.com> * Cubic interpolation fix in arc_01 2D-simulation * Update system tests (#202) * Create new test class, update test_arc01 * Use a config file to organise the tests * Parameterise, and add pytest-check * More tests going in, urls still broken * Update all links to Zenodo sandbox * Remove separate test files - they're all in one now * Update to latest sandbox version * Update to point to published Zenodo * Add docstrings and help * Add pyyaml to requirements.txt Co-authored-by: Sam Cunliffe --- tdms/include/argument_parser.h | 8 + tdms/include/arrays.h | 3 +- tdms/include/cell_coordinate.h | 12 + tdms/include/field.h | 138 +- tdms/include/globals.h | 1 + tdms/include/interpolate.h | 97 -- tdms/include/interpolation_methods.h | 71 +- tdms/include/iterator.h | 3 +- tdms/include/simulation_parameters.h | 6 +- tdms/matlabio/matlabio.h | 3 +- tdms/src/argument_parser.cpp | 5 + tdms/src/fields/base.cpp | 88 +- tdms/src/fields/electric.cpp | 113 +- tdms/src/fields/magnetic.cpp | 263 ++-- tdms/src/interpolate.cpp | 1187 ----------------- tdms/src/interpolation_methods.cpp | 162 +-- tdms/src/iterator.cpp | 300 ++--- tdms/src/main.cpp | 25 +- tdms/src/simulation_parameters.cpp | 4 +- tdms/tests/requirements.txt | 2 + tdms/tests/system/read_config.py | 224 ++++ tdms/tests/system/test_arc01.py | 37 - tdms/tests/system/test_arc02.py | 29 - tdms/tests/system/test_arc03.py | 29 - tdms/tests/system/test_arc08.py | 33 - tdms/tests/system/test_arc09.py | 43 - tdms/tests/system/test_arc10.py | 41 - tdms/tests/system/test_arc12.py | 41 - tdms/tests/system/test_arc13.py | 41 - tdms/tests/system/test_fdtd.py | 49 - tdms/tests/system/test_system.py | 78 ++ tdms/tests/system/utils.py | 31 +- .../field_tests/test_Field_interpolation.cpp | 237 ++-- .../unit/test_interpolation_determination.cpp | 171 ++- .../unit/test_interpolation_functions.cpp | 21 +- 35 files changed, 1207 insertions(+), 2389 deletions(-) delete mode 100644 tdms/include/interpolate.h delete mode 100644 tdms/src/interpolate.cpp create mode 100644 tdms/tests/system/read_config.py delete mode 100644 tdms/tests/system/test_arc01.py delete mode 100644 tdms/tests/system/test_arc02.py delete mode 100644 tdms/tests/system/test_arc03.py delete mode 100644 tdms/tests/system/test_arc08.py delete mode 100644 tdms/tests/system/test_arc09.py delete mode 100644 tdms/tests/system/test_arc10.py delete mode 100644 tdms/tests/system/test_arc12.py delete mode 100644 tdms/tests/system/test_arc13.py delete mode 100644 tdms/tests/system/test_fdtd.py create mode 100644 tdms/tests/system/test_system.py diff --git a/tdms/include/argument_parser.h b/tdms/include/argument_parser.h index a6f0332c4..e6ca32cc1 100644 --- a/tdms/include/argument_parser.h +++ b/tdms/include/argument_parser.h @@ -73,6 +73,14 @@ class ArgumentNamespace{ */ bool finite_difference() const; + /** + * @brief Check whether we were asked to use the cubic interpolation schemes over the BLi schemes (default is no) + * + * @return true if provided '-nbli' or '--no-band-limited' + * @return false otheriwse + */ + bool cubic_interpolation() const; + /** * @brief Gets the input filename * diff --git a/tdms/include/arrays.h b/tdms/include/arrays.h index c9e5353fd..11ec64709 100644 --- a/tdms/include/arrays.h +++ b/tdms/include/arrays.h @@ -10,6 +10,7 @@ #include +#include "cell_coordinate.h" #include "matlabio.h" #include "utils.h" #include "globals.h" @@ -313,7 +314,7 @@ class Tensor3D{ n_rows = _n_rows; } - inline T** operator[] (int value) const { return tensor[value]; }; + inline T **operator[](int value) const { return tensor[value]; }; bool has_elements(){ return tensor != nullptr; }; diff --git a/tdms/include/cell_coordinate.h b/tdms/include/cell_coordinate.h index e2f8d80ef..c86bb0ba0 100644 --- a/tdms/include/cell_coordinate.h +++ b/tdms/include/cell_coordinate.h @@ -3,6 +3,7 @@ * @author William Graham (ccaegra@ucl.ac.uk) * @brief Class declaration for Yee cell coordinates */ +#pragma once class CellCoordinate { @@ -11,11 +12,22 @@ class CellCoordinate { int cell_i = 0, cell_j = 0, cell_k = 0; public: + /** + * @brief Construct a new Cell Coordinate object + * + * @param i,j,k The cell index to assign + */ CellCoordinate(int i = 0, int j = 0, int k = 0) { cell_i = i; cell_j = j; cell_k = k; }; + /** + * @brief Construct a new Cell Coordinate object, passing in a contiguous block of memory + * + * @param indices The indices of the cell index to assign, in the order i,j,k + * @param buffer_start First index to read the indices from + */ CellCoordinate(int *indices, int buffer_start = 0) { cell_i = indices[buffer_start]; cell_j = indices[buffer_start + 1]; diff --git a/tdms/include/field.h b/tdms/include/field.h index 6d9a632b4..b9286a56d 100644 --- a/tdms/include/field.h +++ b/tdms/include/field.h @@ -6,6 +6,7 @@ #include +#include "cell_coordinate.h" #include "dimensions.h" #include "mat_io.h" #include "simulation_parameters.h" @@ -25,18 +26,26 @@ * NOTE: For storage purposes, this means that field values associated to cells are stored _to the left_. * That is, Grid(0,0,0) is associated to the cell (-1,-1,-1). This is contrary to the way values are associated to cells, where cell (0,0,0) is associated to the field values (0,0,0). */ -class Grid{ +class Grid { +protected: + // the preferred interpolation methods (pim) for interpolating between the grid values, default is BandLimited + PreferredInterpolationMethods pim = PreferredInterpolationMethods::BandLimited; public: - int I_tot = 0; - int J_tot = 0; - int K_tot = 0; + int I_tot = 0; + int J_tot = 0; + int K_tot = 0; - /** + /** * Maximum value out of I_tot, J_tot and K_tot * @return value */ - int max_IJK_tot() const {return max(I_tot, J_tot, K_tot); }; + int max_IJK_tot() const { return max(I_tot, J_tot, K_tot); }; + + /** + * @brief Set the preferred interpolation methods + */ + void set_preferred_interpolation_methods(PreferredInterpolationMethods _pim) { pim = _pim; }; }; class SplitFieldComponent: public Tensor3D{ @@ -45,6 +54,9 @@ 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[](CellCoordinate cell) const { return tensor[cell.k()][cell.j()][cell.i()]; } + void initialise_from_matlab(double*** tensor, Dimensions &dims); /** @@ -118,10 +130,10 @@ class SplitField : public Grid{ * @brief Interpolates a SplitField component to the centre of a Yee cell * * @param d SplitField component to interpolate - * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of + * @param cell Index (i,j,k) of the Yee cell to interpolate to the centre of * @return double The interpolated field value */ - virtual double interpolate_to_centre_of(AxialDirection d, int i, int j, int k) = 0; + virtual double interpolate_to_centre_of(AxialDirection d, CellCoordinate cell) = 0; }; class ElectricSplitField: public SplitField{ @@ -142,10 +154,10 @@ class ElectricSplitField: public SplitField{ * @brief Interpolates a split E-field component to the centre of a Yee cell * * @param d Field component to interpolate - * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of + * @param cell Index (i,j,k) of the Yee cell to interpolate to the centre of * @return double The interpolated component value */ - double interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override; + double interpolate_to_centre_of(AxialDirection d, CellCoordinate cell) override; }; class MagneticSplitField: public SplitField{ @@ -166,10 +178,10 @@ class MagneticSplitField: public SplitField{ * @brief Interpolates a split E-field component to the centre of a Yee cell * * @param d Field component to interpolate - * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of + * @param cell Index (i,j,k) of the Yee cell to interpolate to the centre of * @return double The interpolated component value */ - double interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override; + double interpolate_to_centre_of(AxialDirection d, CellCoordinate cell) override; }; class CurrentDensitySplitField: public SplitField{ @@ -186,7 +198,7 @@ class CurrentDensitySplitField: public SplitField{ CurrentDensitySplitField(int I_total, int J_total, int K_total) : SplitField(I_total, J_total, K_total){}; - double interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override { return 0.; }; + double interpolate_to_centre_of(AxialDirection d, CellCoordinate cell) override { return 0.; }; }; /** @@ -194,7 +206,6 @@ class CurrentDensitySplitField: public SplitField{ * at each (x, y, z) grid point */ class Field : public Grid{ - public: double ft = 0.; // TODO: an explanation of what this is @@ -273,15 +284,52 @@ class Field : public Grid{ * @brief Interpolates a Field component to the centre of a Yee cell * * @param d Field component to interpolate - * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of + * @param cell Index (i,j,k) of the Yee cell to interpolate to the centre of * @return std::complex The interpolated field value */ - virtual std::complex interpolate_to_centre_of(AxialDirection d, int i, int j, int k) = 0; + virtual std::complex interpolate_to_centre_of(AxialDirection d, CellCoordinate cell) = 0; + + /** + * @brief Interpolates the Field over the range provided. + * + * Default range is to interpolate to the midpoint of all consecutive points. + * + * @param[out] x_out,y_out,z_out Output arrays for interpolated values + * @param i_lower,j_lower,k_lower Lower index for interpolation in the i,j,k directions, respectively + * @param i_upper,j_upper,k_upper Upper index for interpolation in the i,j,k directions, respectively + * @param mode Determines which field components to compute, based on the simulation Dimension + */ + void interpolate_over_range(mxArray **x_out, mxArray **y_out, mxArray **z_out, int i_lower, + int i_upper, int j_lower, int j_upper, int k_lower, int k_upper, + Dimension mode = Dimension::THREE); + void interpolate_over_range(mxArray **x_out, mxArray **y_out, mxArray **z_out, + Dimension mode = Dimension::THREE); + + /** + * @brief Interpolates the Field's transverse electric components to the centre of Yee cell i,j,k + * + * @param[in] cell Yee cell index + * @param[out] x_at_centre,y_at_centre,z_at_centre Addresses to write interpolated values for the x,y,z components (respectively) + */ + virtual void interpolate_transverse_electric_components(CellCoordinate cell, + std::complex *x_at_centre, + std::complex *y_at_centre, + std::complex *z_at_centre) = 0; + /** + * @brief Interpolates the Field's transverse magnetic components to the centre of Yee cell i,j,k + * + * @param[in] cell Yee cell index + * @param[out] x_at_centre,y_at_centre,z_at_centre Addresses to write interpolated values for the x,y,z components (respectively) + */ + virtual void interpolate_transverse_magnetic_components(CellCoordinate cell, + std::complex *x_at_centre, + std::complex *y_at_centre, + std::complex *z_at_centre) = 0; /** * Set the values of all components in this field from another, equally sized field */ - void set_values_from(Field &other); + void set_values_from(Field &other); ~Field(); }; @@ -302,7 +350,32 @@ class ElectricField: public Field{ * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of * @return std::complex The interpolated component value */ - std::complex interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override; + std::complex interpolate_to_centre_of(AxialDirection d, CellCoordinate cell) override; + + /** + * @brief Interpolates the transverse electric components to the centre of Yee cell i,j,k. + * + * Ex and Ey are interpolated. Ez is set to a placeholder (default) value. + * + * @param[in] cell Yee cell index + * @param[out] x_at_centre,y_at_centre,z_at_centre Addresses to write interpolated values for the x,y,z components (respectively) + */ + void interpolate_transverse_electric_components(CellCoordinate cell, + std::complex *x_at_centre, + std::complex *y_at_centre, + std::complex *z_at_centre) override; + /** + * @brief Interpolates the transverse magnetic components to the centre of Yee cell i,j,k. + * + * Ez is interpolated. Ex and Ey are set to a placeholder (default) values. + * + * @param[in] cell Yee cell index + * @param[out] x_at_centre,y_at_centre,z_at_centre Addresses to write interpolated values for the x,y,z components (respectively) + */ + void interpolate_transverse_magnetic_components(CellCoordinate cell, + std::complex *x_at_centre, + std::complex *y_at_centre, + std::complex *z_at_centre) override; }; class MagneticField: public Field{ @@ -318,10 +391,35 @@ class MagneticField: public Field{ * @brief Interpolates an H-field component to the centre of a Yee cell * * @param d Field component to interpolate - * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of + * @param cell Index (i,j,k) of the Yee cell to interpolate to the centre of * @return std::complex The interpolated component value */ - std::complex interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override; + std::complex interpolate_to_centre_of(AxialDirection d, CellCoordinate cell) override; + + /** + * @brief Interpolates the transverse electric components to the centre of Yee cell i,j,k. + * + * Hz is interpolated. Hx and Hy are set to a placeholder (default) values. + * + * @param[in] cell Yee cell index + * @param[out] x_at_centre,y_at_centre,z_at_centre Addresses to write interpolated values for the x,y,z components (respectively) + */ + void interpolate_transverse_electric_components(CellCoordinate cell, + std::complex *x_at_centre, + std::complex *y_at_centre, + std::complex *z_at_centre) override; + /** + * @brief Interpolates the transverse magnetic components to the centre of Yee cell i,j,k. + * + * Hx and Hy are interpolated. Hz is set to a placeholder (default) value. + * + * @param[in] cell Yee cell index + * @param[out] x_at_centre,y_at_centre,z_at_centre Addresses to write interpolated values for the x,y,z components (respectively) + */ + void interpolate_transverse_magnetic_components(CellCoordinate cell, + std::complex *x_at_centre, + std::complex *y_at_centre, + std::complex *z_at_centre) override; }; /** diff --git a/tdms/include/globals.h b/tdms/include/globals.h index 5cf1365f9..0267c4a03 100644 --- a/tdms/include/globals.h +++ b/tdms/include/globals.h @@ -53,6 +53,7 @@ enum AxialDirection enum ModeOfRun { Pass1 , Pass2 }; enum RCSType { parallel , perpendicular }; enum SolverMethod { PseudoSpectral, FiniteDifference }; +enum PreferredInterpolationMethods { BandLimited, Cubic }; // ************************************** // Mathematical Constants diff --git a/tdms/include/interpolate.h b/tdms/include/interpolate.h deleted file mode 100644 index 969c8b08f..000000000 --- a/tdms/include/interpolate.h +++ /dev/null @@ -1,97 +0,0 @@ -/** - * @file interpolate.h - * @brief Interpolation of field values within FDTD grid. - * - * The Yee cell specifies the 6 field components at different points in space. - * Interpolation is required to calculate the 6 field components at the same - * spatial position. - */ -#pragma once - -#include - -#include "field.h" - -/** - * @brief Interpolate the electric field to the origin of the Yee cell - * - * @param[in] Ex_yee Steady state x component of electric field calculated at points in the Yee cell - * @param[in] Ey_yee Steady state y component of electric field calculated at points in the Yee cell - * @param[in] Ez_yee Steady state z component of electric field calculated at points in the Yee cell - * @param[out] Ex Steady state x component of electric field interpolated to Yee cell origin - * @param[out] Ey Steady state y component of electric field interpolated to Yee cell origin - * @param[out] Ez Steady state z component of electric field interpolated to Yee cell origin - * @param[in] I Number of elements in the i direction of the FDTD grid - * @param[in] J Number of elements in the j direction of the FDTD grid - * @param[in] K Number of elements in the k direction of the FDTD grid - * @param[in] i_l Least i index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest j index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest k index into the FDTD grid to evaluate the field at. Should be <= K-2 - */ -void interpolateFieldCentralE( double ***Ex_yee, double ***Ey_yee, double ***Ez_yee, - double ***Ex , double ***Ey , double ***Ez , - int I , int J , int K , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u); -void interpolateFieldCentralE_TE( double ***Ex_yee, double ***Ey_yee, double ***Ez_yee, - double ***Ex , double ***Ey , double ***Ez , - int I , int J , int K , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u); -void mxInterpolateFieldCentralE( mxArray *Ex_yee , mxArray *Ey_yee , mxArray *Ez_yee, - mxArray **Ex , mxArray **Ey , mxArray **Ez , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u); -void mxInterpolateFieldCentralE_TE( mxArray *Ex_yee , mxArray *Ey_yee , mxArray *Ez_yee, - mxArray **Ex , mxArray **Ey , mxArray **Ez , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u); -void mxInterpolateFieldCentralE_TM( mxArray *Ex_yee , mxArray *Ey_yee , mxArray *Ez_yee, - mxArray **Ex , mxArray **Ey , mxArray **Ez , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u); -void interpolateFieldCentralH( double ***Hx_yee, double ***Hy_yee, double ***Hz_yee, - double ***Hx , double ***Hy , double ***Hz , - int I , int J , int K , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u); -void interpolateFieldCentralH_TE( double ***Hx_yee, double ***Hy_yee, double ***Hz_yee, - double ***Hx , double ***Hy , double ***Hz , - int I , int J , int K , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u); -void mxInterpolateFieldCentralH( mxArray *Hx_yee , mxArray *Hy_yee , mxArray *Hz_yee, - mxArray **Hx , mxArray **Hy , mxArray **Hz , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u); -void mxInterpolateFieldCentralH_TE( mxArray *Hx_yee , mxArray *Hy_yee , mxArray *Hz_yee, - mxArray **Hx , mxArray **Hy , mxArray **Hz , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u); -void mxInterpolateFieldCentralH_TM( mxArray *Hx_yee , mxArray *Hy_yee , mxArray *Hz_yee, - mxArray **Hx , mxArray **Hy , mxArray **Hz , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u); -void interpolateTimeDomainFieldCentralE( SplitFieldComponent& Exy, SplitFieldComponent& Exz, SplitFieldComponent& Eyx, SplitFieldComponent& Eyz, SplitFieldComponent& Ezx, SplitFieldComponent& Ezy, - int i, int j, int k, - double *Ex, double *Ey, double *Ez); -void interpolateTimeDomainFieldCentralEBandLimited( SplitFieldComponent& Exy, SplitFieldComponent& Exz, SplitFieldComponent& Eyx, SplitFieldComponent& Eyz, SplitFieldComponent& Ezx, SplitFieldComponent& Ezy, - int i, int j, int k, - double *Ex, double *Ey, double *Ez); -void interpolateTimeDomainFieldCentralE_2Dy( SplitFieldComponent& Exy, SplitFieldComponent& Exz, SplitFieldComponent& Eyx, SplitFieldComponent& Eyz, SplitFieldComponent& Ezx, SplitFieldComponent& Ezy, - int i, int j, int k, - double *Ex, double *Ey, double *Ez); -void interpolateTimeDomainFieldCentralE_TE( SplitFieldComponent& Exy, SplitFieldComponent& Exz, SplitFieldComponent& Eyx, SplitFieldComponent& Eyz, SplitFieldComponent& Ezx, SplitFieldComponent& Ezy, - int i, int j, int k, - double *Ex, double *Ey, double *Ez); -void interpolateTimeDomainFieldCentralE_TM( SplitFieldComponent& Exy, SplitFieldComponent& Exz, SplitFieldComponent& Eyx, SplitFieldComponent& Eyz, SplitFieldComponent& Ezx, SplitFieldComponent& Ezy, - int i, int j, int k, - double *Ex, double *Ey, double *Ez); -void interpolateTimeDomainFieldCentralH( SplitFieldComponent& Hxy, SplitFieldComponent& Hxz, SplitFieldComponent& Hyx, SplitFieldComponent& Hyz, SplitFieldComponent& Hzx, SplitFieldComponent& Hzy, - int i, int j, int k, - double *Hx, double *Hy, double *Hz); -void interpolateTimeDomainFieldCentralH_2Dy( SplitFieldComponent& Hxy, SplitFieldComponent& Hxz, SplitFieldComponent& Hyx, SplitFieldComponent& Hyz, SplitFieldComponent& Hzx, SplitFieldComponent& Hzy, - int i, int j, int k, - double *Hx, double *Hy, double *Hz); -void interpolateTimeDomainFieldCentralH_TE( SplitFieldComponent& Hxy, SplitFieldComponent& Hxz, SplitFieldComponent& Hyx, SplitFieldComponent& Hyz, SplitFieldComponent& Hzx, SplitFieldComponent& Hzy, - int i, int j, int k, - double *Hx, double *Hy, double *Hz); -void interpolateTimeDomainFieldCentralH_TM( SplitFieldComponent& Hxy, SplitFieldComponent& Hxz, SplitFieldComponent& Hyx, SplitFieldComponent& Hyz, SplitFieldComponent& Hzx, SplitFieldComponent& Hzy, - int i, int j, int k, - double *Hx, double *Hy, double *Hz); -void interpolateTimeDomainFieldCentralHBandLimited( SplitFieldComponent& Hxy, SplitFieldComponent& Hxz, SplitFieldComponent& Hyx, SplitFieldComponent& Hyz, SplitFieldComponent& Hzx, SplitFieldComponent& Hzy, - int i, int j, int k, - double *Hx, double *Hy, double *Hz); diff --git a/tdms/include/interpolation_methods.h b/tdms/include/interpolation_methods.h index 0086efeb1..a13a0ecba 100644 --- a/tdms/include/interpolation_methods.h +++ b/tdms/include/interpolation_methods.h @@ -9,53 +9,7 @@ #include -/*Use cubic interpolation to interpolate between the middle 2 of 4 points - * v0 v1 v2 v3 - * o o x o o - */ -double interp1(double v1, double v2, double v3, double v4); -/*Use cubic interpolation to interpolate between the middle 2 of 4 points - * v0 v1 v2 v3 - * o o x o o - */ -double interp1(double *v); -/*Use cubic interpolation to interpolate between the first 2 of 4 points - * v0 v1 v2 v3 - * o x o o o - */ -double interp2(double v1, double v2, double v3, double v4); -/*Use cubic interpolation to interpolate between the first 2 of 4 points - * v0 v1 v2 v3 - * o x o o o - */ -double interp2(double *v); -/*Use cubic interpolation to interpolate between the last 2 of 4 points - * v0 v1 v2 v3 - * o o o x o - */ -double interp3(double v1, double v2, double v3, double v4); -/*Use cubic interpolation to interpolate between the last 2 of 4 points - * v0 v1 v2 v3 - * o o o x o - */ -double interp3(double *v); - -/** - * @brief Checks whether the limits of field extraction are within range of the FDTD grid - * - * Since cubic interpolation is being used, it must be ensured that (in any given direction) that the least index used is no smaller than 2, whilst the greatest no larger than the maximum number of cells in that direction - 2. - * - * @param i_l Least i index into the FDTD grid to evaluate the field at - * @param i_u Greatest i index into the FDTD grid to evaluate the field at - * @param j_l Least j index into the FDTD grid to evaluate the field at - * @param j_u Greatest j index into the FDTD grid to evaluate the field at - * @param k_l Least k index into the FDTD grid to evaluate the field at - * @param k_u Greatest k index into the FDTD grid to evaluate the field at - * @param nI,nJ,nK Number of elements in the i,j,k (respectively) direction of the FDTD grid - * - * @throws runtime_error In the event that the limits of field extraction are outside the FDTD grid - */ -void checkInterpolationPoints(int i_l, int i_u, int j_l, int j_u, int k_l, int k_u, int nI, int nJ, int nK); +#include "globals.h" /** * @brief Defines our order of preference for the use of the various schemes. @@ -79,11 +33,11 @@ enum scheme_value BAND_LIMITED_4 = 7, // use bandlimited interpolation w/ interp position = 4 BAND_LIMITED_5 = 5, // use bandlimited interpolation w/ interp position = 5 BAND_LIMITED_6 = 3, // use bandlimited interpolation w/ interp position = 6 - BAND_LIMITED_7 = -1, // use bandlimited interpolation w/ interp position = 7 [Only applicable if we want to extend beyond the final Yee cell, current code functionality is to throw an error in the case where this would be used.] - BAND_LIMITED_CELL_ZERO = -2, // use bandlimited interpolation to interpolate to the centre of Yee cell <0 [implemented, but current code functionality is to throw an error here] - CUBIC_INTERP_MIDDLE = 2, // cubic interpolation to middle 2 of 4 points (interp1) - CUBIC_INTERP_FIRST = 1, // cubic interpolation to first 2 of 4 points (interp2) - CUBIC_INTERP_LAST = 0 // cubic interpolation to last 2 of 4 points (interp3) + BAND_LIMITED_7 = 2, // use bandlimited interpolation w/ interp position = 7 + BAND_LIMITED_CELL_ZERO = 1, // use bandlimited interpolation to interpolate to the centre of Yee cell 0 + CUBIC_INTERP_MIDDLE = 0, // cubic interpolation to middle 2 of 4 points (interp1) + CUBIC_INTERP_FIRST = -1, // cubic interpolation to first 2 of 4 points (interp2) + CUBIC_INTERP_LAST = -2 // cubic interpolation to last 2 of 4 points (interp3) }; class InterpolationScheme { @@ -175,10 +129,13 @@ const InterpolationScheme CBMid = InterpolationScheme(CUBIC_INTERP_MIDDLE); const InterpolationScheme CBLst = InterpolationScheme(CUBIC_INTERP_LAST); /** - * @brief Determines the appropriate interpolation scheme to use, given the current cell and number of cells in a given dimension. + * @brief Determine the appropriate interpolation scheme to use, given the number of datapoints available and the position to which we wish to interpolate to. * - * @param cells_in_direction The number of cells in the direction parallel to interpolation - * @param cell_id The current cell index to interpolate to the centre of - * @return const interpScheme* The interpolation scheme that should be used + * @param datapts_in_direction The number of datapoints available along this axis + * @param interpolation_position Interpolation is to be performed to the midpoint of the datapoints indexed with (interpolation_position-1) and (interpolation_position) + * @param interpolation_methods The preferred interpolation methods to use + * @return const InterpolationScheme& */ -const InterpolationScheme &best_scheme(int cells_in_direction, int cell_id); +const InterpolationScheme &best_scheme(int datapts_in_direction, int interpolation_position, + PreferredInterpolationMethods interpolation_methods = + PreferredInterpolationMethods::BandLimited); diff --git a/tdms/include/iterator.h b/tdms/include/iterator.h index 3176c9d8e..561c900b5 100644 --- a/tdms/include/iterator.h +++ b/tdms/include/iterator.h @@ -17,7 +17,8 @@ * @brief Executes the main simulation. * Used to be the MATLAB mexFunction */ -void execute_simulation(int, mxArray**, int, const mxArray**, SolverMethod); +void execute_simulation(int, mxArray **, int, const mxArray **, SolverMethod solver_method, + PreferredInterpolationMethods preferred_interpolation_methods); void extractPhasorsPlane( double **iwave_lEx_Rbs, double **iwave_lEx_Ibs, double **iwave_lEy_Rbs, double **iwave_lEy_Ibs, double **iwave_lHx_Rbs, double **iwave_lHx_Ibs, double **iwave_lHy_Rbs, double **iwave_lHy_Ibs, diff --git a/tdms/include/simulation_parameters.h b/tdms/include/simulation_parameters.h index 5c7d6a77d..8f811a28d 100644 --- a/tdms/include/simulation_parameters.h +++ b/tdms/include/simulation_parameters.h @@ -34,9 +34,9 @@ enum RunMode{ }; enum Dimension{ - THREE, //< Full dimensionality - compute all H and E components - TE, //< Transverse electric - only compute Ex, Ey, and Hz components - TM //< Transverse magnetic - only compute Hx, Hy, and Ez components + THREE, //< Full dimensionality - compute all H and E components + TRANSVERSE_ELECTRIC, //< Transverse electric - only compute Ex, Ey, and Hz components + TRANSVERSE_MAGNETIC //< Transverse magnetic - only compute Hx, Hy, and Ez components }; enum InterpolationMethod{ diff --git a/tdms/matlabio/matlabio.h b/tdms/matlabio/matlabio.h index f24814138..53d48e35e 100644 --- a/tdms/matlabio/matlabio.h +++ b/tdms/matlabio/matlabio.h @@ -31,7 +31,7 @@ void free_cast_matlab_4D_array(double ****castArray, int nlayers, int nblocks); /** * Casts a 3-dimensional array such that it may be indexed according to the usual array indexing * scheme array[k,j,i]. - * @param array is a pointer to a matlab 4 dimensional array + * @param array is a pointer to a matlab 3 dimensional array * @param nrows the number of rows in the array * @param ncols the number of columns in the array @@ -39,7 +39,6 @@ void free_cast_matlab_4D_array(double ****castArray, int nlayers, int nblocks); */ template T*** cast_matlab_3D_array(T *array, S nrows, S ncols, S nlayers){ - T ***p; nlayers = std::max(nlayers, S(1)); p = (T ***)malloc((unsigned) (nlayers*sizeof(T **))); diff --git a/tdms/src/argument_parser.cpp b/tdms/src/argument_parser.cpp index a95e8c724..1e6c096e0 100644 --- a/tdms/src/argument_parser.cpp +++ b/tdms/src/argument_parser.cpp @@ -52,6 +52,7 @@ void ArgumentParser::print_help_message(){ "Options:\n" "-h:\tDisplay this help message\n" "-fd, --finite-difference:\tUse the finite-difference solver, instead of the pseudo-spectral method.\n" + "-c, --cubic-interpolation:\tUse cubic interpolation to determine field values at Yee cell centres, as opposed to band-limited interpolation.\n" "-q:\tQuiet operation. Silence all logging\n" "-m:\tMinimise output file size by not saving vertex and facet information\n\n"); } @@ -87,6 +88,10 @@ bool ArgumentNamespace::finite_difference() const { return this->have_flag("-fd") || this->have_flag("--finite-difference"); } +bool ArgumentNamespace::cubic_interpolation() const { + return this->have_flag("-c") || this->have_flag("--cubic-interpolation"); +} + const char* ArgumentNamespace::output_filename() { if (has_grid_filename()){ diff --git a/tdms/src/fields/base.cpp b/tdms/src/fields/base.cpp index 5c50dd189..cd4403806 100644 --- a/tdms/src/fields/base.cpp +++ b/tdms/src/fields/base.cpp @@ -10,8 +10,6 @@ void Field::add_to_angular_norm(int n, int Nt, SimulationParameters ¶ms) { angular_norm += phasor_norm(ft, n, params.omega_an, params.dt, Nt); } - - void Field::normalise_volume() { double norm_r = std::real(angular_norm); @@ -38,7 +36,7 @@ Field::Field(int I_total, int J_total, int K_total) { void Field::allocate() { for (auto arr : {&real, &imag}) { - arr->allocate(I_tot+1, J_tot+1, K_tot+1); + arr->allocate(I_tot, J_tot, K_tot); } }; @@ -57,7 +55,91 @@ complex Field::phasor_norm(double f, int n, double omega, double dt, int * 1./((double) Nt); } +void Field::interpolate_over_range(mxArray **x_out, mxArray **y_out, mxArray **z_out, int i_lower, + int i_upper, int j_lower, int j_upper, int k_lower, int k_upper, + Dimension mode) { + // prepare output dimensions + const int ndims = 3; + int outdims[ndims] = {i_upper - i_lower + 1, j_upper - j_lower + 1, 1}; + if (mode == THREE) { + outdims[2] = k_upper - k_lower + 1; + if (outdims[1] < 1) { + // full simulation (all components being computed) but in 2D - allow one cell in the y-direction to avoid NULL pointers + outdims[1] = 1; + } + } + // create the memory for the outputs + *x_out = mxCreateNumericArray(ndims, (const mwSize *) outdims, mxDOUBLE_CLASS, mxCOMPLEX); + *y_out = mxCreateNumericArray(ndims, (const mwSize *) outdims, mxDOUBLE_CLASS, mxCOMPLEX); + *z_out = mxCreateNumericArray(ndims, (const mwSize *) outdims, mxDOUBLE_CLASS, mxCOMPLEX); + // allow memory to cast to our datatypes + XYZTensor3D real_out, imag_out; + real_out.x = cast_matlab_3D_array(mxGetPr(*x_out), outdims[0], outdims[1], outdims[2]); + imag_out.x = cast_matlab_3D_array(mxGetPi(*x_out), outdims[0], outdims[1], outdims[2]); + real_out.y = cast_matlab_3D_array(mxGetPr(*y_out), outdims[0], outdims[1], outdims[2]); + imag_out.y = cast_matlab_3D_array(mxGetPi(*y_out), outdims[0], outdims[1], outdims[2]); + real_out.z = cast_matlab_3D_array(mxGetPr(*z_out), outdims[0], outdims[1], outdims[2]); + imag_out.z = cast_matlab_3D_array(mxGetPi(*z_out), outdims[0], outdims[1], outdims[2]); + + /* If we are not interpolating _all_ components, we need to use E- and H-field specific methods. + Also (in full-field 3D simulations only) we need to check whether j_upper x_at_centre = interpolate_to_centre_of(AxialDirection::X, current_cell), + y_at_centre = interpolate_to_centre_of(AxialDirection::Y, current_cell), + z_at_centre = interpolate_to_centre_of(AxialDirection::Z, current_cell); + real_out.x[k - k_lower][0][i - i_lower] = x_at_centre.real(); + imag_out.x[k - k_lower][0][i - i_lower] = x_at_centre.imag(); + // y interpolation doesn't take place, so use placeholder values + real_out.y[k - k_lower][0][i - i_lower] = y_at_centre.real(); + imag_out.y[k - k_lower][0][i - i_lower] = y_at_centre.imag(); + real_out.z[k - k_lower][0][i - i_lower] = z_at_centre.real(); + imag_out.z[k - k_lower][0][i - i_lower] = z_at_centre.imag(); + } + } + } else { + // 3D loop, and use the field-specific methods to ensure we interpolate components correctly + for (int i = i_lower; i <= i_upper; i++) { + for (int j = j_lower; j <= j_upper; j++) { + for (int k = k_lower; k <= k_upper; k++) { + complex x_at_centre, y_at_centre, z_at_centre; + CellCoordinate current_cell(i,j,k); + switch (mode) { + case Dimension::THREE: + // 3D interpolation is identical for both E and H fields + x_at_centre = interpolate_to_centre_of(AxialDirection::X, current_cell); + y_at_centre = interpolate_to_centre_of(AxialDirection::Y, current_cell); + z_at_centre = interpolate_to_centre_of(AxialDirection::Z, current_cell); + break; + case Dimension::TRANSVERSE_ELECTRIC: + interpolate_transverse_electric_components(current_cell, &x_at_centre, &y_at_centre, &z_at_centre); + break; + case Dimension::TRANSVERSE_MAGNETIC: + interpolate_transverse_magnetic_components(current_cell, &x_at_centre, &y_at_centre, &z_at_centre); + break; + } + real_out.x[k - k_lower][j - j_lower][i - i_lower] = x_at_centre.real(); + imag_out.x[k - k_lower][j - j_lower][i - i_lower] = x_at_centre.imag(); + real_out.y[k - k_lower][j - j_lower][i - i_lower] = y_at_centre.real(); + imag_out.y[k - k_lower][j - j_lower][i - i_lower] = y_at_centre.imag(); + real_out.z[k - k_lower][j - j_lower][i - i_lower] = z_at_centre.real(); + imag_out.z[k - k_lower][j - j_lower][i - i_lower] = z_at_centre.imag(); + } + } + } + } +} +void Field::interpolate_over_range(mxArray **x_out, mxArray **y_out, mxArray **z_out, + Dimension mode) { + interpolate_over_range(x_out, y_out, z_out, 0, I_tot, 0, J_tot, 0, K_tot, mode); +} Field::~Field() { diff --git a/tdms/src/fields/electric.cpp b/tdms/src/fields/electric.cpp index 845b46b9b..643fa5409 100644 --- a/tdms/src/fields/electric.cpp +++ b/tdms/src/fields/electric.cpp @@ -10,49 +10,66 @@ double ElectricField::phase(int n, double omega, double dt){ return omega * ((double) n + 1) * dt; } -/* -You may notice several +1's appearing when the interpolation data (interp_data) arrays are constructed from the field values. -This is because E.x[k][j][i] is, by construction of the Grid class, not the value associated to Yee cell (i,j,k) but rather cell (i-1,j,k). Similarly Ey[k][j][i] is the value of cell (i,j-1,k) and Ez[k][j][i] of cell (i,j,k-1). -*/ +void ElectricField::interpolate_transverse_electric_components(CellCoordinate cell, complex *x_at_centre, + complex *y_at_centre, + complex *z_at_centre) { + *x_at_centre = interpolate_to_centre_of(AxialDirection::X, cell); + *y_at_centre = interpolate_to_centre_of(AxialDirection::Y, cell); + *z_at_centre = complex(0., 0.); +} +void ElectricField::interpolate_transverse_magnetic_components(CellCoordinate cell, complex *x_at_centre, + complex *y_at_centre, + complex *z_at_centre) { + *x_at_centre = complex(0., 0.); + *y_at_centre = complex(0., 0.); + *z_at_centre = interpolate_to_centre_of(AxialDirection::Z, cell); +} -complex ElectricField::interpolate_to_centre_of(AxialDirection d, int i, int j, int k) { +complex ElectricField::interpolate_to_centre_of(AxialDirection d, CellCoordinate cell) { const InterpolationScheme *scheme; + int i = cell.i(), j = cell.j(), k = cell.k(); // prepare input data - if using a cubic scheme we have reserved more memory than necessary but nevermind complex interp_data[8]; + switch (d) { case X: // determine the interpolation scheme to use - scheme = &(best_scheme(I_tot, i)); + scheme = &(best_scheme(I_tot, i, pim)); // now fill the interpolation data // i - (scheme.number_of_datapoints_to_left) is the index of the Yee cell that plays the role of v0 in the interpolation for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { interp_data[ind] = - real.x[k][j][i + 1 - scheme->number_of_datapoints_to_left + ind] + - IMAGINARY_UNIT * imag.x[k][j][i + 1 - scheme->number_of_datapoints_to_left + ind]; + real.x[k][j][i - scheme->number_of_datapoints_to_left + ind] + + IMAGINARY_UNIT * imag.x[k][j][i - scheme->number_of_datapoints_to_left + ind]; } break; case Y: - // determine the interpolation scheme to use - scheme = &(best_scheme(J_tot, j)); + // 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 (J_tot <= 1) { + return complex(real.y[k][0][i], imag.y[k][0][i]); + } else { // 3D simulation, interpolation is as normal + // determine the interpolation scheme to use + scheme = &(best_scheme(J_tot, j, pim)); - // now fill the interpolation data - // j - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation - for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { - interp_data[ind] = - real.y[k][j + 1 - scheme->number_of_datapoints_to_left + ind][i] + - IMAGINARY_UNIT * imag.y[k][j + 1 - scheme->number_of_datapoints_to_left + ind][i]; + // now fill the interpolation data + // j - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation + for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { + interp_data[ind] = + real.y[k][j - scheme->number_of_datapoints_to_left + ind][i] + + IMAGINARY_UNIT * imag.y[k][j - scheme->number_of_datapoints_to_left + ind][i]; + } } break; case Z: // determine the interpolation scheme to use - scheme = &(best_scheme(K_tot, k)); + scheme = &(best_scheme(K_tot, k, pim)); // now fill the interpolation data // k - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { interp_data[ind] = - real.z[k + 1 - scheme->number_of_datapoints_to_left + ind][j][i] + - IMAGINARY_UNIT * imag.z[k + 1 - scheme->number_of_datapoints_to_left + ind][j][i]; + real.z[k - scheme->number_of_datapoints_to_left + ind][j][i] + + IMAGINARY_UNIT * imag.z[k - scheme->number_of_datapoints_to_left + ind][j][i]; } break; default: @@ -63,47 +80,53 @@ complex ElectricField::interpolate_to_centre_of(AxialDirection d, int i, return scheme->interpolate(interp_data); } -double ElectricSplitField::interpolate_to_centre_of(AxialDirection d, int i, int j, int k) { +double ElectricSplitField::interpolate_to_centre_of(AxialDirection d, CellCoordinate cell) { const InterpolationScheme *scheme; + int i = cell.i(), j = cell.j(), k = cell.k(); // prepare input data - if using a cubic scheme we have reserved more memory than necessary but nevermind double interp_data[8]; switch (d) { case X: - scheme = &(best_scheme(I_tot, i)); + scheme = &(best_scheme(I_tot, i, pim)); // now fill the interpolation data // i - (scheme.number_of_datapoints_to_left) is the index of the Yee cell that plays the role of v0 in the interpolation for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { - interp_data[ind] = xy[k][j][i + 1 - scheme->number_of_datapoints_to_left + ind] + - xz[k][j][i + 1 - scheme->number_of_datapoints_to_left + 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]; } // now run the interpolation scheme and place the result into the output return scheme->interpolate(interp_data); break; case Y: - scheme = &(best_scheme(J_tot, j)); - // now fill the interpolation data - // j - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation - for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { - interp_data[ind] = yx[k][j + 1 - scheme->number_of_datapoints_to_left + ind][i] + - yz[k][j + 1 - scheme->number_of_datapoints_to_left + ind][i]; + // 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 (J_tot <= 1) { + return yx[k][0][i] + yz[k][0][i]; + } else {// 3D simulation, interpolation is as normal + scheme = &(best_scheme(J_tot, j, pim)); + // now fill the interpolation data + // j - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation + 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]; + } + // now run the interpolation scheme and place the result into the output + return scheme->interpolate(interp_data); + break; + case Z: + scheme = &(best_scheme(K_tot, k, pim)); + // now fill the interpolation data + // k - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation + 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]; + } + // now run the interpolation scheme and place the result into the output + return scheme->interpolate(interp_data); } - // now run the interpolation scheme and place the result into the output - return scheme->interpolate(interp_data); break; - case Z: - scheme = &(best_scheme(K_tot, k)); - // now fill the interpolation data - // k - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation - for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { - interp_data[ind] = zx[k + 1 - scheme->number_of_datapoints_to_left + ind][j][i] + - zy[k + 1 - scheme->number_of_datapoints_to_left + ind][j][i]; + default: + throw runtime_error("Invalid axial direction selected for interpolation!\n"); + break; } - // now run the interpolation scheme and place the result into the output - return scheme->interpolate(interp_data); - break; - default: - throw runtime_error("Invalid axial direction selected for interpolation!\n"); - break; - } } diff --git a/tdms/src/fields/magnetic.cpp b/tdms/src/fields/magnetic.cpp index b59215e00..1f6535e68 100644 --- a/tdms/src/fields/magnetic.cpp +++ b/tdms/src/fields/magnetic.cpp @@ -1,5 +1,6 @@ #include "field.h" +#include "cell_coordinate.h" #include "interpolation_methods.h" #include "globals.h" @@ -10,6 +11,23 @@ double MagneticField::phase(int n, double omega, double dt){ return omega * ((double) n + 0.5) * dt; // 0.5 added because it's known half a time step after E } +void MagneticField::interpolate_transverse_electric_components(CellCoordinate cell, + complex *x_at_centre, + complex *y_at_centre, + complex *z_at_centre) { + *x_at_centre = complex(0., 0.); + *y_at_centre = complex(0., 0.); + *z_at_centre = interpolate_to_centre_of(AxialDirection::Z, cell); +} +void MagneticField::interpolate_transverse_magnetic_components(CellCoordinate cell, + complex *x_at_centre, + complex *y_at_centre, + complex *z_at_centre) { + *x_at_centre = interpolate_to_centre_of(AxialDirection::X, cell); + *y_at_centre = interpolate_to_centre_of(AxialDirection::Y, cell); + *z_at_centre = complex(0., 0.); +} + /* 2D INTERPOLATION SCHEMES (FOR THE MAGNETIC FIELD IN 3D SIMULATIONS) Unlike the E-field, the H-field components associated with Yee cell i,j,k are _not_ aligned with the centre of the Yee cells. @@ -46,79 +64,89 @@ b_j - b_scheme.number_of_datapoints_to_left + b_scheme.first_nonzero_coeff <= ce Now with approximations of Ha at each (a_i, cell_b + Db, c_k), we can pass this information to b_scheme.interpolate to recover the value of Ha at the centre of Yee cell (a_i, b_j, c_k). In a 2D simulation, J_tot = 0. This means we cannot interpolate in two directions for the Hx and Hz fields, since there is no y-direction to interpolate in. As a result, we only interpolate in the z-direction for Hx and x-direction for Hz when running a two-dimensional simulation. - -You may notice several +1's appearing when the interpolation data (interp_data) arrays are constructed from the field values. -This is because H.x[k][j][i] is, by construction of the Grid class, not the value associated to Yee cell (i,j,k) but rather cell (i,j-1,k-1). Similarly Ey[k][j][i] is the value of cell (i-1,j,k) and Ez[k][j][i] of cell (i-1,j-1,k). */ -complex MagneticField::interpolate_to_centre_of(AxialDirection d, int i, int j, int k) { +complex MagneticField::interpolate_to_centre_of(AxialDirection d, CellCoordinate cell) { // this data will be passed to the second interpolation scheme complex data_for_second_scheme[8]; // this data will hold values for the interpolation in the first interpolation scheme complex data_for_first_scheme[8]; // the interpolation schemes that are to be used const InterpolationScheme *b_scheme, *c_scheme; + int i = cell.i(), j = cell.j(), k = cell.k(); switch (d) { case X: // Associations: a = x, b = y, c = z - c_scheme = &(best_scheme(K_tot, k)); - b_scheme = &(best_scheme(J_tot, j)); - - if (c_scheme->is_better_than(*b_scheme)) { - // we will be interpolating in the z-direction first, then in y - for (int jj = b_scheme->first_nonzero_coeff; jj <= b_scheme->last_nonzero_coeff; jj++) { - // this is the j-index of the cell we are looking at - int cell_j = j + 1 - b_scheme->number_of_datapoints_to_left + jj; - // determine the Hx values of cells (i, cell_j, k-z_scheme.index-1) through (i, cell_j, k-z_scheme.index-1+7), and interpolate them via z_scheme - for (int kk = c_scheme->first_nonzero_coeff; kk <= c_scheme->last_nonzero_coeff; kk++) { - // the k-index of the current cell we are looking at (readability, don't need to define this here) - int cell_k = k + 1 - c_scheme->number_of_datapoints_to_left + kk; - // gather the data for interpolating in the z dimension - data_for_first_scheme[kk] = - real.x[cell_k][cell_j][i] + IMAGINARY_UNIT * imag.x[cell_k][cell_j][i]; - } - // 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 data being passed to the y_scheme - data_for_second_scheme[jj] = c_scheme->interpolate(data_for_first_scheme); + if (J_tot <=1) { + // this is a 2D simulation + // we simply interpolate the Hx field in the z-direction to get the field value at the centre + c_scheme = &(best_scheme(K_tot, k, pim)); + for(int kk = c_scheme->first_nonzero_coeff; kk <= c_scheme->last_nonzero_coeff; kk++) { + int cell_k = k - c_scheme->number_of_datapoints_to_left + kk; + data_for_first_scheme[kk] = real.x[cell_k][0][i] + IMAGINARY_UNIT * imag.x[cell_k][0][i]; } - // now interpolate in the y-direction to the centre of Yee cell (i,j,k) - return b_scheme->interpolate(data_for_second_scheme); - } else { - // we will be interpolating in the y-direction first, then in z - for (int kk = c_scheme->first_nonzero_coeff; kk <= c_scheme->last_nonzero_coeff; kk++) { - // this is the k-index of the cell we are looking at - int cell_k = k + 1 - c_scheme->number_of_datapoints_to_left + kk; - // determine the Hx values of cells (i, j - y_scheme.index-1, cell_k) through (i, j - y_scheme.index-1+7, cell_k), and interpolate them via y_scheme + return c_scheme->interpolate(data_for_first_scheme); + } + else { + c_scheme = &(best_scheme(K_tot, k, pim)); + b_scheme = &(best_scheme(J_tot, j, pim)); + + if (c_scheme->is_better_than(*b_scheme)) { + // we will be interpolating in the z-direction first, then in y for (int jj = b_scheme->first_nonzero_coeff; jj <= b_scheme->last_nonzero_coeff; jj++) { - // the j-index of the current cell we are looking at (readability, don't need to define this here) - int cell_j = j + 1 - b_scheme->number_of_datapoints_to_left + jj; - // gather the data for interpolating in the y dimension - data_for_first_scheme[jj] = - real.x[cell_k][cell_j][i] + IMAGINARY_UNIT * imag.x[cell_k][cell_j][i]; + // this is the j-index of the cell we are looking at + int cell_j = j - b_scheme->number_of_datapoints_to_left + jj; + // determine the Hx values of cells (i, cell_j, k-z_scheme.index-1) through (i, cell_j, k-z_scheme.index-1+7), and interpolate them via z_scheme + for (int kk = c_scheme->first_nonzero_coeff; kk <= c_scheme->last_nonzero_coeff; kk++) { + // the k-index of the current cell we are looking at (readability, don't need to define this here) + 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] = + real.x[cell_k][cell_j][i] + IMAGINARY_UNIT * imag.x[cell_k][cell_j][i]; + } + // 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 data being passed to the y_scheme + data_for_second_scheme[jj] = c_scheme->interpolate(data_for_first_scheme); } - // 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 data being passed to the y_scheme - data_for_second_scheme[kk] = b_scheme->interpolate(data_for_first_scheme); + // now interpolate in the y-direction to the centre of Yee cell (i,j,k) + return b_scheme->interpolate(data_for_second_scheme); + } else { + // we will be interpolating in the y-direction first, then in z + for (int kk = c_scheme->first_nonzero_coeff; kk <= c_scheme->last_nonzero_coeff; kk++) { + // this is the k-index of the cell we are looking at + int cell_k = k - c_scheme->number_of_datapoints_to_left + kk; + // determine the Hx values of cells (i, j - y_scheme.index-1, cell_k) through (i, j - y_scheme.index-1+7, cell_k), and interpolate them via y_scheme + for (int jj = b_scheme->first_nonzero_coeff; jj <= b_scheme->last_nonzero_coeff; jj++) { + // the j-index of the current cell we are looking at (readability, don't need to define this here) + 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] = + real.x[cell_k][cell_j][i] + IMAGINARY_UNIT * imag.x[cell_k][cell_j][i]; + } + // 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 data being passed to the y_scheme + data_for_second_scheme[kk] = b_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the z-direction to the centre of Yee cell (i,j,k) + return c_scheme->interpolate(data_for_second_scheme); } - // now interpolate in the z-direction to the centre of Yee cell (i,j,k) - return c_scheme->interpolate(data_for_second_scheme); } break; case Y: // Associations: a = y, b = z, c = x - c_scheme = &(best_scheme(I_tot, i)); - b_scheme = &(best_scheme(K_tot, k)); + c_scheme = &(best_scheme(I_tot, i, pim)); + b_scheme = &(best_scheme(K_tot, k, pim)); if (b_scheme->is_better_than(*c_scheme)) { // we will be interpolating in the z-direction first, then in x for (int ii = c_scheme->first_nonzero_coeff; ii <= c_scheme->last_nonzero_coeff; ii++) { // this is the i-index of the cell we are looking at - int cell_i = i + 1 - c_scheme->number_of_datapoints_to_left + ii; + int cell_i = i - c_scheme->number_of_datapoints_to_left + ii; // determine the Hy values of cells (cell_i, j, k-z_scheme.index-1) through (cell_i, j, k-z_scheme.index-1+7), and interpolate them via z_scheme for (int kk = b_scheme->first_nonzero_coeff; kk <= b_scheme->last_nonzero_coeff; kk++) { // the k-index of the current cell we are looking at (readability, don't need to define this here) - int cell_k = k + 1 - b_scheme->number_of_datapoints_to_left + kk; + 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] = real.y[cell_k][j][cell_i] + IMAGINARY_UNIT * imag.y[cell_k][j][cell_i]; @@ -133,11 +161,11 @@ complex MagneticField::interpolate_to_centre_of(AxialDirection d, int i, // we will be interpolating in the x-direction first, then in z for (int kk = b_scheme->first_nonzero_coeff; kk <= b_scheme->last_nonzero_coeff; kk++) { // this is the k-index of the cell we are looking at - int cell_k = k + 1 - b_scheme->number_of_datapoints_to_left + kk; + int cell_k = k - b_scheme->number_of_datapoints_to_left + kk; // determine the Hy values of cells (i - x_scheme.index-1, j, cell_k) through (i- x_scheme.index-1+7, j, cell_k), and interpolate them via x_scheme for (int ii = c_scheme->first_nonzero_coeff; ii <= c_scheme->last_nonzero_coeff; ii++) { // the i-index of the current cell we are looking at (readability, don't need to define this here) - int cell_i = i + 1 - c_scheme->number_of_datapoints_to_left + ii; + 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] = real.y[cell_k][j][cell_i] + IMAGINARY_UNIT * imag.y[cell_k][j][cell_i]; @@ -152,47 +180,61 @@ complex MagneticField::interpolate_to_centre_of(AxialDirection d, int i, break; case Z: // Associations: a = z, b = x, c = y - b_scheme = &(best_scheme(I_tot, i)); - c_scheme = &(best_scheme(J_tot, j)); - - if (c_scheme->is_better_than(*b_scheme)) { - // we will be interpolating in the y-direction first, then in x + if (J_tot <= 1) { + // this is a 2D simulation + // we simply interpolate the Hx field in the z-direction to get the field value at the centre + b_scheme = &(best_scheme(I_tot, i, pim)); for (int ii = b_scheme->first_nonzero_coeff; ii <= b_scheme->last_nonzero_coeff; ii++) { - // this is the i-index of the cell we are looking at - int cell_i = i + 1 - b_scheme->number_of_datapoints_to_left + ii; - // determine the Hz values of cells (cell_i, j-y_scheme.index-1, k) through (cell_i, j-y_scheme.index-1+7, k), and interpolate them via y_scheme - for (int jj = c_scheme->first_nonzero_coeff; jj <= c_scheme->last_nonzero_coeff; jj++) { - // the j-index of the current cell we are looking at (readability, don't need to define this here) - int cell_j = j + 1 - c_scheme->number_of_datapoints_to_left + jj; - // gather the data for interpolating in the y dimension - data_for_first_scheme[jj] = - real.z[k][cell_j][cell_i] + IMAGINARY_UNIT * imag.z[k][cell_j][cell_i]; - } - // 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 data being passed to the x_scheme - data_for_second_scheme[ii] = c_scheme->interpolate(data_for_first_scheme); + int cell_i = i - b_scheme->number_of_datapoints_to_left + ii; + data_for_first_scheme[ii] = real.z[k][0][cell_i] + IMAGINARY_UNIT * imag.z[k][0][cell_i]; } - // now interpolate in the x-direction to the centre of Yee cell (i,j,k) - return b_scheme->interpolate(data_for_second_scheme); - } else { - // we will be interpolating in the x-direction first, then in y - for (int jj = c_scheme->first_nonzero_coeff; jj <= c_scheme->last_nonzero_coeff; jj++) { - // this is the j-index of the cell we are looking at - int cell_j = j + 1 - c_scheme->number_of_datapoints_to_left + jj; - // determine the Hz values of cells (i - x_scheme.index-1, cell_j, k) through (i- x_scheme.index-1+7, cell_j, k), and interpolate them via x_scheme + return b_scheme->interpolate(data_for_first_scheme); + } + else { + b_scheme = &(best_scheme(I_tot, i, pim)); + c_scheme = &(best_scheme(J_tot, j, pim)); + + if (c_scheme->is_better_than(*b_scheme)) { + // we will be interpolating in the y-direction first, then in x for (int ii = b_scheme->first_nonzero_coeff; ii <= b_scheme->last_nonzero_coeff; ii++) { - // the i-index of the current cell we are looking at (readability, don't need to define this here) - int cell_i = i + 1 - b_scheme->number_of_datapoints_to_left + ii; - // gather the data for interpolating in the x dimension - data_for_first_scheme[ii] = - real.z[k][cell_j][cell_i] + IMAGINARY_UNIT * imag.z[k][cell_j][cell_i]; + // this is the i-index of the cell we are looking at + int cell_i = i - b_scheme->number_of_datapoints_to_left + ii; + // determine the Hz values of cells (cell_i, j-y_scheme.index-1, k) through (cell_i, j-y_scheme.index-1+7, k), and interpolate them via y_scheme + for (int jj = c_scheme->first_nonzero_coeff; jj <= c_scheme->last_nonzero_coeff; jj++) { + // the j-index of the current cell we are looking at (readability, don't need to define this here) + 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] = + real.z[k][cell_j][cell_i] + IMAGINARY_UNIT * imag.z[k][cell_j][cell_i]; + } + + // 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 data being passed to the x_scheme + data_for_second_scheme[ii] = c_scheme->interpolate(data_for_first_scheme); } - // 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 data being passed to the y_scheme - data_for_second_scheme[jj] = b_scheme->interpolate(data_for_first_scheme); + // now interpolate in the x-direction to the centre of Yee cell (i,j,k) + return b_scheme->interpolate(data_for_second_scheme); + } else { + // we will be interpolating in the x-direction first, then in y + for (int jj = c_scheme->first_nonzero_coeff; jj <= c_scheme->last_nonzero_coeff; jj++) { + // this is the j-index of the cell we are looking at + int cell_j = j - c_scheme->number_of_datapoints_to_left + jj; + // determine the Hz values of cells (i - x_scheme.index-1, cell_j, k) through (i- x_scheme.index-1+7, cell_j, k), and interpolate them via x_scheme + for (int ii = b_scheme->first_nonzero_coeff; ii <= b_scheme->last_nonzero_coeff; ii++) { + // the i-index of the current cell we are looking at (readability, don't need to define this here) + 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] = + real.z[k][cell_j][cell_i] + IMAGINARY_UNIT * imag.z[k][cell_j][cell_i]; + } + + // 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 data being passed to the y_scheme + data_for_second_scheme[jj] = b_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the y-direction to the centre of Yee cell (i,j,k) + return c_scheme->interpolate(data_for_second_scheme); } - // now interpolate in the y-direction to the centre of Yee cell (i,j,k) - return c_scheme->interpolate(data_for_second_scheme); } break; default: @@ -201,8 +243,9 @@ complex MagneticField::interpolate_to_centre_of(AxialDirection d, int i, } } -double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int i, int j, int k) { +double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, CellCoordinate cell) { const InterpolationScheme *b_scheme, *c_scheme; + int i = cell.i(), j = cell.j(), k = cell.k(); // this data will be passed to the second interpolation scheme double data_for_second_scheme[8]; // this data will hold values for the interpolation in the first interpolation scheme @@ -210,33 +253,33 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int i, int switch (d) { case X: - if (J_tot <= 0) { + if (J_tot <= 1) { // in a 2D simulation, we must interpolate in the z-direction to recover Hx, due to the magnetic-field offsets from the centre - b_scheme = &(best_scheme(K_tot, k)); + b_scheme = &(best_scheme(K_tot, k, pim)); // now fill the interpolation data // i - (scheme.number_of_datapoints_to_left) is the index of the Yee cell that plays the role of v0 in the interpolation for (int ind = b_scheme->first_nonzero_coeff; ind <= b_scheme->last_nonzero_coeff; ind++) { data_for_first_scheme[ind] = - xy[k + 1 - b_scheme->number_of_datapoints_to_left + ind][j][i] + - xz[k + 1 - b_scheme->number_of_datapoints_to_left + ind][j][i]; + xy[k - b_scheme->number_of_datapoints_to_left + ind][j][i] + + xz[k - b_scheme->number_of_datapoints_to_left + ind][j][i]; } // now run the interpolation scheme and place the result into the output return b_scheme->interpolate(data_for_first_scheme); } else { // Associations: a = x, b = y, c = z - c_scheme = &(best_scheme(K_tot, k)); - b_scheme = &(best_scheme(J_tot, j)); + c_scheme = &(best_scheme(K_tot, k, pim)); + b_scheme = &(best_scheme(J_tot, j, pim)); if (c_scheme->is_better_than(*b_scheme)) { // we will be interpolating in the z-direction first, then in y for (int jj = b_scheme->first_nonzero_coeff; jj <= b_scheme->last_nonzero_coeff; jj++) { // this is the j-index of the cell we are looking at - int cell_j = j + 1 - b_scheme->number_of_datapoints_to_left + jj; + int cell_j = j - b_scheme->number_of_datapoints_to_left + jj; // determine the Hx values of cells (i, cell_j, k-z_scheme.index-1) through (i, cell_j, k-z_scheme.index-1+7), and interpolate them via z_scheme for (int kk = c_scheme->first_nonzero_coeff; kk <= c_scheme->last_nonzero_coeff; kk++) { // the k-index of the current cell we are looking at (readability, don't need to define this here) - int cell_k = k + 1 - c_scheme->number_of_datapoints_to_left + kk; + 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]; } @@ -250,11 +293,11 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int i, int // we will be interpolating in the y-direction first, then in z for (int kk = c_scheme->first_nonzero_coeff; kk <= c_scheme->last_nonzero_coeff; kk++) { // this is the k-index of the cell we are looking at - int cell_k = k + 1 - c_scheme->number_of_datapoints_to_left + kk; + int cell_k = k - c_scheme->number_of_datapoints_to_left + kk; // determine the Hx values of cells (i, j - y_scheme.index-1, cell_k) through (i, j - y_scheme.index-1+7, cell_k), and interpolate them via y_scheme for (int jj = b_scheme->first_nonzero_coeff; jj <= b_scheme->last_nonzero_coeff; jj++) { // the j-index of the current cell we are looking at (readability, don't need to define this here) - int cell_j = j + 1 - b_scheme->number_of_datapoints_to_left + jj; + 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]; } @@ -270,18 +313,18 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int i, int case Y: // we can always interpolate in two directions for Hy, since even in a 2D simulation the x- and z-directions still exist // Associations: a = y, b = z, c = x - c_scheme = &(best_scheme(I_tot, i)); - b_scheme = &(best_scheme(K_tot, k)); + c_scheme = &(best_scheme(I_tot, i, pim)); + b_scheme = &(best_scheme(K_tot, k, pim)); if (b_scheme->is_better_than(*c_scheme)) { // we will be interpolating in the z-direction first, then in x for (int ii = c_scheme->first_nonzero_coeff; ii <= c_scheme->last_nonzero_coeff; ii++) { // this is the i-index of the cell we are looking at - int cell_i = i + 1 - c_scheme->number_of_datapoints_to_left + ii; + int cell_i = i - c_scheme->number_of_datapoints_to_left + ii; // determine the Hy values of cells (cell_i, j, k-z_scheme.index-1) through (cell_i, j, k-z_scheme.index-1+7), and interpolate them via z_scheme for (int kk = b_scheme->first_nonzero_coeff; kk <= b_scheme->last_nonzero_coeff; kk++) { // the k-index of the current cell we are looking at (readability, don't need to define this here) - int cell_k = k + 1 - b_scheme->number_of_datapoints_to_left + kk; + 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]; } @@ -295,11 +338,11 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int i, int // we will be interpolating in the x-direction first, then in z for (int kk = b_scheme->first_nonzero_coeff; kk <= b_scheme->last_nonzero_coeff; kk++) { // this is the k-index of the cell we are looking at - int cell_k = k + 1 - b_scheme->number_of_datapoints_to_left + kk; + int cell_k = k - b_scheme->number_of_datapoints_to_left + kk; // determine the Hy values of cells (i - x_scheme.index-1, j, cell_k) through (i- x_scheme.index-1+7, j, cell_k), and interpolate them via x_scheme for (int ii = c_scheme->first_nonzero_coeff; ii <= c_scheme->last_nonzero_coeff; ii++) { // the i-index of the current cell we are looking at (readability, don't need to define this here) - int cell_i = i + 1 - c_scheme->number_of_datapoints_to_left + ii; + 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]; } @@ -312,33 +355,33 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int i, int } break; case Z: - if (J_tot <= 0) { + if (J_tot <= 1) { // in a 2D simulation, we must interpolate in the x-direction to recover Hz, due to the magnetic-field offsets from the centre - b_scheme = &(best_scheme(I_tot, i)); + b_scheme = &(best_scheme(I_tot, i, pim)); // now fill the interpolation data // i - (scheme.number_of_datapoints_to_left) is the index of the Yee cell that plays the role of v0 in the interpolation for (int ind = b_scheme->first_nonzero_coeff; ind <= b_scheme->last_nonzero_coeff; ind++) { data_for_first_scheme[ind] = - zx[k][j][i + 1 - b_scheme->number_of_datapoints_to_left + ind] + - zy[k][j][i + 1 - b_scheme->number_of_datapoints_to_left + 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]; } // now run the interpolation scheme and place the result into the output return b_scheme->interpolate(data_for_first_scheme); } else { // Associations: a = z, b = x, c = y - b_scheme = &(best_scheme(I_tot, i)); - c_scheme = &(best_scheme(J_tot, j)); + b_scheme = &(best_scheme(I_tot, i, pim)); + c_scheme = &(best_scheme(J_tot, j, pim)); if (c_scheme->is_better_than(*b_scheme)) { // we will be interpolating in the y-direction first, then in x for (int ii = b_scheme->first_nonzero_coeff; ii <= b_scheme->last_nonzero_coeff; ii++) { // this is the i-index of the cell we are looking at - int cell_i = i + 1 - b_scheme->number_of_datapoints_to_left + ii; + int cell_i = i - b_scheme->number_of_datapoints_to_left + ii; // determine the Hz values of cells (cell_i, j-y_scheme.index-1, k) through (cell_i, j-y_scheme.index-1+7, k), and interpolate them via y_scheme for (int jj = c_scheme->first_nonzero_coeff; jj <= c_scheme->last_nonzero_coeff; jj++) { // the j-index of the current cell we are looking at (readability, don't need to define this here) - int cell_j = j + 1 - c_scheme->number_of_datapoints_to_left + jj; + 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]; } @@ -352,11 +395,11 @@ double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int i, int // we will be interpolating in the x-direction first, then in y for (int jj = c_scheme->first_nonzero_coeff; jj <= c_scheme->last_nonzero_coeff; jj++) { // this is the j-index of the cell we are looking at - int cell_j = j + 1 - c_scheme->number_of_datapoints_to_left + jj; + int cell_j = j - c_scheme->number_of_datapoints_to_left + jj; // determine the Hz values of cells (i - x_scheme.index-1, cell_j, k) through (i- x_scheme.index-1+7, cell_j, k), and interpolate them via x_scheme for (int ii = b_scheme->first_nonzero_coeff; ii <= b_scheme->last_nonzero_coeff; ii++) { // the i-index of the current cell we are looking at (readability, don't need to define this here) - int cell_i = i + 1 - b_scheme->number_of_datapoints_to_left + ii; + 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]; } diff --git a/tdms/src/interpolate.cpp b/tdms/src/interpolate.cpp deleted file mode 100644 index 4c290dac1..000000000 --- a/tdms/src/interpolate.cpp +++ /dev/null @@ -1,1187 +0,0 @@ -#include "interpolate.h" - -#include - -#include "interpolation_methods.h" -#include "matlabio.h" - -using namespace std; - - -void interpolateFieldCentralE( double ***Ex_yee, double ***Ey_yee, double ***Ez_yee, - double ***Ex , double ***Ey , double ***Ez , - int I , int J , int K , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - int i,j,k; - - //first check that the limits of field extraction of within range - checkInterpolationPoints(i_l, i_u, j_l, j_u, k_l, k_u, I, J, K); - - //fprintf(stderr, "interpolateFieldCentralE 01\n"); - if(j_u= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest j index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest k index into the FDTD grid to evaluate the field at. Should be <= K-2 - */ -void interpolateFieldCentralE_TE( double ***Ex_yee, double ***Ey_yee, double ***Ez_yee, - double ***Ex , double ***Ey , double ***Ez , - int I , int J , int K , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - int i,j,k; - - //first check that the limits of field extraction of within range - checkInterpolationPoints(i_l, i_u, j_l, j_u, k_l, k_u, I, J, K); - - for(k=k_l;k<=k_u;k++) - for(j=j_l;j<=j_u;j++) - for(i=i_l;i<=i_u;i++){ - Ex[k-k_l][j-j_l][i-i_l] = interp1(Ex_yee[k][j][i-2], Ex_yee[k][j][i-1], Ex_yee[k][j][i], Ex_yee[k][j][i+1]); - Ey[k-k_l][j-j_l][i-i_l] = interp1(Ey_yee[k][j-2][i], Ey_yee[k][j-1][i], Ey_yee[k][j][i], Ey_yee[k][j+1][i]); - Ez[k-k_l][j-j_l][i-i_l] = 0.; - } -} - -/** - * @brief Interpolate the TM electric field (z component) to the origin of the Yee cell - * - * @param[in] Ex_yee Steady state x component of electric field calculated at points in the Yee cell - * @param[in] Ey_yee Steady state y component of electric field calculated at points in the Yee cell - * @param[in] Ez_yee Steady state z component of electric field calculated at points in the Yee cell - * @param[out] Ex Steady state x component of electric field interpolated to Yee cell origin - * @param[out] Ey Steady state y component of electric field interpolated to Yee cell origin - * @param[out] Ez Steady state z component of electric field interpolated to Yee cell origin - * @param[in] I Number of elements in the i direction of the FDTD grid - * @param[in] J Number of elements in the j direction of the FDTD grid - * @param[in] K Number of elements in the k direction of the FDTD grid - * @param[in] i_l Least i index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest j index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest k index into the FDTD grid to evaluate the field at. Should be <= K-2 - */ -void interpolateFieldCentralE_TM( double ***Ex_yee, double ***Ey_yee, double ***Ez_yee, - double ***Ex , double ***Ey , double ***Ez , - int I , int J , int K , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - int i,j,k; - - //first check that the limits of field extraction of within range - checkInterpolationPoints(i_l, i_u, j_l, j_u, k_l, k_u, I, J, K); - - for(k=k_l;k<=k_u;k++) - for(j=j_l;j<=j_u;j++) - for(i=i_l;i<=i_u;i++){ - Ex[k-k_l][j-j_l][i-i_l] = 0.; - Ey[k-k_l][j-j_l][i-i_l] = 0.; - Ez[k-k_l][j-j_l][i-i_l] = Ez_yee[k][j][i]; - } -} - -/** - * @brief Interpolate the electric field to the origin of the Yee cell - MATLAB interface - * - * @param[in] Ex_yee Input Yee cell Ex field matrix - * @param[in] Ey_yee Input Yee cell Ey field matrix - * @param[in] Ez_yee Input Yee cell Ez field matrix - * @param[out] Ex Output Yee cell Ex field matrix - * @param[out] Ey Output Yee cell Ey field matrix - * @param[out] Ez Output Yee cell Ez field matrix - * @param[in] i_l Least i index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= K-2 - * - * I, J, K are the number of elements in the i, j, k directions of the FDTD grid, respectively. - */ -void mxInterpolateFieldCentralE( mxArray *Ex_yee , mxArray *Ey_yee , mxArray *Ez_yee, - mxArray **Ex , mxArray **Ey , mxArray **Ez , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - double ***ExR, ***ExI, ***EyR, ***EyI, ***EzR, ***EzI,***Ex_yee_R, ***Ex_yee_I, ***Ey_yee_R, ***Ey_yee_I, ***Ez_yee_R, ***Ez_yee_I; - const int *indims; - int outdims[3]; - - //fprintf(stderr, "mxInterpolateFieldCentralE Pos 00\n"); - if( (int)mxGetNumberOfDimensions( (const mxArray *)Ex_yee) < 3){ - throw runtime_error("Error in mxInterpolateFieldCentralE, Ex_yee does not have 3 dimensions\n"); - } - //fprintf(stderr, "mxInterpolateFieldCentralE Pos 01\n"); - - indims = (int *)mxGetDimensions( (mxArray *)Ex_yee); - //fprintf(stderr, "mxInterpolateFieldCentralE(indims): (%d,%d,%d)\n",indims[0],indims[1],indims[2]); - //fprintf(stderr, "mxInterpolateFieldCentralE: j_u: %d, j_l: %d\n",j_u,j_l); - //assume that all matrices have the same dimensions - DANGEROUS, SHOULDN'T WE CHECK THIS?? - if( !mxIsComplex(Ex_yee) ){ - mexErrMsgTxt("Ex_yee is not complex"); - } - if( !mxIsComplex(Ey_yee) ){ - mexErrMsgTxt("Ey_yee is not complex"); - } - if( !mxIsComplex(Ez_yee) ){ - mexErrMsgTxt("Ez_yee is not complex"); - } - - // cast data in the Yee cells to complex MATLAB arrays - //fprintf(stderr, "mxInterpolateFieldCentralE Pos 02\n"); - Ex_yee_R = cast_matlab_3D_array(mxGetPr(Ex_yee), indims[0], indims[1], indims[2]); - Ex_yee_I = cast_matlab_3D_array(mxGetPi(Ex_yee), indims[0], indims[1], indims[2]); - - Ey_yee_R = cast_matlab_3D_array(mxGetPr(Ey_yee), indims[0], indims[1], indims[2]); - Ey_yee_I = cast_matlab_3D_array(mxGetPi(Ey_yee), indims[0], indims[1], indims[2]); - - Ez_yee_R = cast_matlab_3D_array(mxGetPr(Ez_yee), indims[0], indims[1], indims[2]); - Ez_yee_I = cast_matlab_3D_array(mxGetPi(Ez_yee), indims[0], indims[1], indims[2]); - //fprintf(stderr, "mxInterpolateFieldCentralE Pos 03\n"); - - //now construct the output matrices - int ndims = 3; - - outdims[0] = i_u - i_l + 1; - outdims[1] = j_u - j_l + 1; - if(outdims[1]<1) - outdims[1]=1; - outdims[2] = k_u - k_l + 1; - //fprintf(stderr, "mxInterpolateFieldCentralE(outdims): (%d,%d,%d)\n",outdims[0],outdims[1],outdims[2]); - //fprintf(stderr, "mxInterpolateFieldCentralE Pos 04\n"); - *Ex = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Ey = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Ez = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); -//fprintf(stderr, "mxInterpolateFieldCentralE Pos 04a\n"); - ExR = cast_matlab_3D_array(mxGetPr(*Ex), outdims[0], outdims[1], outdims[2]); - ExI = cast_matlab_3D_array(mxGetPi(*Ex), outdims[0], outdims[1], outdims[2]); -//fprintf(stderr, "mxInterpolateFieldCentralE Pos 04b\n"); - EyR = cast_matlab_3D_array(mxGetPr(*Ey), outdims[0], outdims[1], outdims[2]); - EyI = cast_matlab_3D_array(mxGetPi(*Ey), outdims[0], outdims[1], outdims[2]); -//fprintf(stderr, "mxInterpolateFieldCentralE Pos 04c\n"); - EzR = cast_matlab_3D_array(mxGetPr(*Ez), outdims[0], outdims[1], outdims[2]); - EzI = cast_matlab_3D_array(mxGetPi(*Ez), outdims[0], outdims[1], outdims[2]); - //fprintf(stderr, "mxInterpolateFieldCentralE Pos 05\n"); - - //now interpolate fields - // imaginary part - interpolateFieldCentralE( Ex_yee_I, Ey_yee_I, Ez_yee_I, ExI, EyI, EzI, indims[0], indims[1], indims[2], i_l, i_u, j_l, j_u, k_l, k_u); - // real part - interpolateFieldCentralE( Ex_yee_R, Ey_yee_R, Ez_yee_R, ExR, EyR, EzR, indims[0], indims[1], indims[2], i_l, i_u, j_l, j_u, k_l, k_u); - - //free the extra memory used by casting array - free_cast_matlab_3D_array(Ex_yee_R, indims[2]); - free_cast_matlab_3D_array(Ex_yee_I, indims[2]); - free_cast_matlab_3D_array(Ey_yee_R, indims[2]); - free_cast_matlab_3D_array(Ey_yee_I, indims[2]); - free_cast_matlab_3D_array(Ez_yee_R, indims[2]); - free_cast_matlab_3D_array(Ez_yee_I, indims[2]); - - free_cast_matlab_3D_array(ExR, outdims[2]); - free_cast_matlab_3D_array(ExI, outdims[2]); - free_cast_matlab_3D_array(EyR, outdims[2]); - free_cast_matlab_3D_array(EyI, outdims[2]); - free_cast_matlab_3D_array(EzR, outdims[2]); - free_cast_matlab_3D_array(EzI, outdims[2]); -} - -/** - * @brief Interpolate the TE electric field (x,y components) to the origin of the Yee cell - MATLAB interface - * - * @param[in] Ex_yee Input Yee cell Ex field matrix - * @param[in] Ey_yee Input Yee cell Ey field matrix - * @param[in] Ez_yee Input Yee cell Ez field matrix - * @param[out] Ex Output Yee cell Ex field matrix - * @param[out] Ey Output Yee cell Ey field matrix - * @param[out] Ez Output Yee cell Ez field matrix - * @param[in] i_l Least i index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= K-2 - * - * I, J, K are the number of elements in the i, j, k directions of the FDTD grid, respectively. - */ -void mxInterpolateFieldCentralE_TE( mxArray *Ex_yee , mxArray *Ey_yee , mxArray *Ez_yee, - mxArray **Ex , mxArray **Ey , mxArray **Ez , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - double ***ExR, ***ExI, ***EyR, ***EyI, ***EzR, ***EzI,***Ex_yee_R, ***Ex_yee_I, ***Ey_yee_R, ***Ey_yee_I, ***Ez_yee_R, ***Ez_yee_I; - const int *indims; - int outdims[3], ndims; - - if( mxGetNumberOfDimensions(Ex_yee) !=2 ){ - throw runtime_error("Error in mxInterpolateFieldCentralE_TE, Ex_yee does not have 2 dimensions\n"); - - - } - - indims = (int *)mxGetDimensions( (mxArray *)Ex_yee); - //assume that all matrices have the same dimensions - if( !mxIsComplex(Ex_yee) ){ - mexErrMsgTxt("Ex_yee is not complex"); - } - if( !mxIsComplex(Ey_yee) ){ - mexErrMsgTxt("Ey_yee is not complex"); - } - if( !mxIsComplex(Ez_yee) ){ - mexErrMsgTxt("Ez_yee is not complex"); - } - - Ex_yee_R = cast_matlab_3D_array(mxGetPr(Ex_yee), indims[0], indims[1], 0); - Ex_yee_I = cast_matlab_3D_array(mxGetPi(Ex_yee), indims[0], indims[1], 0); - - Ey_yee_R = cast_matlab_3D_array(mxGetPr(Ey_yee), indims[0], indims[1], 0); - Ey_yee_I = cast_matlab_3D_array(mxGetPi(Ey_yee), indims[0], indims[1], 0); - - Ez_yee_R = cast_matlab_3D_array(mxGetPr(Ez_yee), indims[0], indims[1], 0); - Ez_yee_I = cast_matlab_3D_array(mxGetPi(Ez_yee), indims[0], indims[1], 0); - - //now construct the output matrices - ndims = 3; - - outdims[0] = i_u - i_l + 1; - outdims[1] = j_u - j_l + 1; - outdims[2] = 1; - - *Ex = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Ey = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Ez = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - - ExR = cast_matlab_3D_array(mxGetPr(*Ex), outdims[0], outdims[1], outdims[2]); - ExI = cast_matlab_3D_array(mxGetPi(*Ex), outdims[0], outdims[1], outdims[2]); - - EyR = cast_matlab_3D_array(mxGetPr(*Ey), outdims[0], outdims[1], outdims[2]); - EyI = cast_matlab_3D_array(mxGetPi(*Ey), outdims[0], outdims[1], outdims[2]); - - EzR = cast_matlab_3D_array(mxGetPr(*Ez), outdims[0], outdims[1], outdims[2]); - EzI = cast_matlab_3D_array(mxGetPi(*Ez), outdims[0], outdims[1], outdims[2]); - //now finally ready for interpolation - interpolateFieldCentralE_TE( Ex_yee_I, Ey_yee_I, Ez_yee_I, - ExI , EyI , EzI , - indims[0] , indims[1] , indims[2] , - i_l, i_u, j_l, j_u, k_l, k_u); - - interpolateFieldCentralE_TE( Ex_yee_R, Ey_yee_R, Ez_yee_R, - ExR , EyR , EzR , - indims[0] , indims[1] , indims[2] , - i_l, i_u, j_l, j_u, k_l, k_u); - //free the extra memory used by casting array - free_cast_matlab_3D_array(Ex_yee_R, 0); - free_cast_matlab_3D_array(Ex_yee_I, 0); - free_cast_matlab_3D_array(Ey_yee_R, 0); - free_cast_matlab_3D_array(Ey_yee_I, 0); - free_cast_matlab_3D_array(Ez_yee_R, 0); - free_cast_matlab_3D_array(Ez_yee_I, 0); - - free_cast_matlab_3D_array(ExR, outdims[2]); - free_cast_matlab_3D_array(ExI, outdims[2]); - free_cast_matlab_3D_array(EyR, outdims[2]); - free_cast_matlab_3D_array(EyI, outdims[2]); - free_cast_matlab_3D_array(EzR, outdims[2]); - free_cast_matlab_3D_array(EzI, outdims[2]); -} - -/** - * @brief Interpolate the TM electric field (z component) to the origin of the Yee cell - MATLAB interface - * - * @param[in] Ex_yee Input Yee cell Ex field matrix - * @param[in] Ey_yee Input Yee cell Ey field matrix - * @param[in] Ez_yee Input Yee cell Ez field matrix - * @param[out] Ex Output Yee cell Ex field matrix - * @param[out] Ey Output Yee cell Ey field matrix - * @param[out] Ez Output Yee cell Ez field matrix - * @param[in] i_l Least i index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= K-2 - * - * I, J, K are the number of elements in the i, j, k directions of the FDTD grid, respectively. - */ -void mxInterpolateFieldCentralE_TM( mxArray *Ex_yee , mxArray *Ey_yee , mxArray *Ez_yee, - mxArray **Ex , mxArray **Ey , mxArray **Ez , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - double ***ExR, ***ExI, ***EyR, ***EyI, ***EzR, ***EzI,***Ex_yee_R, ***Ex_yee_I, ***Ey_yee_R, ***Ey_yee_I, ***Ez_yee_R, ***Ez_yee_I; - const int *indims; - int outdims[3], ndims; - - if( mxGetNumberOfDimensions(Ex_yee) != 2){ - throw runtime_error("Error in mxInterpolateFieldCentralE_TM, Ex_yee does not have 2 dimensions\n"); - - - } - - indims = (int *)mxGetDimensions(Ex_yee); - //assume that all matrices have the same dimensions - if( !mxIsComplex(Ex_yee) ){ - mexErrMsgTxt("Ex_yee is not complex"); - } - if( !mxIsComplex(Ey_yee) ){ - mexErrMsgTxt("Ey_yee is not complex"); - } - if( !mxIsComplex(Ez_yee) ){ - mexErrMsgTxt("Ez_yee is not complex"); - } - - - Ex_yee_R = cast_matlab_3D_array(mxGetPr(Ex_yee), indims[0], indims[1], 0); - Ex_yee_I = cast_matlab_3D_array(mxGetPi(Ex_yee), indims[0], indims[1], 0); - - Ey_yee_R = cast_matlab_3D_array(mxGetPr(Ey_yee), indims[0], indims[1], 0); - Ey_yee_I = cast_matlab_3D_array(mxGetPi(Ey_yee), indims[0], indims[1], 0); - - Ez_yee_R = cast_matlab_3D_array(mxGetPr(Ez_yee), indims[0], indims[1], 0); - Ez_yee_I = cast_matlab_3D_array(mxGetPi(Ez_yee), indims[0], indims[1], 0); - - //now construct the output matrices - ndims = 3; - - outdims[0] = i_u - i_l + 1; - outdims[1] = j_u - j_l + 1; - outdims[2] = 1; - - *Ex = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Ey = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Ez = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - - ExR = cast_matlab_3D_array(mxGetPr(*Ex), outdims[0], outdims[1], outdims[2]); - ExI = cast_matlab_3D_array(mxGetPi(*Ex), outdims[0], outdims[1], outdims[2]); - - EyR = cast_matlab_3D_array(mxGetPr(*Ey), outdims[0], outdims[1], outdims[2]); - EyI = cast_matlab_3D_array(mxGetPi(*Ey), outdims[0], outdims[1], outdims[2]); - - EzR = cast_matlab_3D_array(mxGetPr(*Ez), outdims[0], outdims[1], outdims[2]); - EzI = cast_matlab_3D_array(mxGetPi(*Ez), outdims[0], outdims[1], outdims[2]); - - //now finally ready for interpolation - interpolateFieldCentralE_TM( Ex_yee_I, Ey_yee_I, Ez_yee_I, - ExI , EyI , EzI , - indims[0] , indims[1] , indims[2] , - i_l, i_u, j_l, j_u, k_l, k_u); - - - interpolateFieldCentralE_TM( Ex_yee_R, Ey_yee_R, Ez_yee_R, - ExR , EyR , EzR , - indims[0] , indims[1] , indims[2] , - i_l, i_u, j_l, j_u, k_l, k_u); - - //free the extra memory used by casting array - free_cast_matlab_3D_array(Ex_yee_R, 0); - free_cast_matlab_3D_array(Ex_yee_I, 0); - free_cast_matlab_3D_array(Ey_yee_R, 0); - free_cast_matlab_3D_array(Ey_yee_I, 0); - free_cast_matlab_3D_array(Ez_yee_R, 0); - free_cast_matlab_3D_array(Ez_yee_I, 0); - - free_cast_matlab_3D_array(ExR, outdims[2]); - free_cast_matlab_3D_array(ExI, outdims[2]); - free_cast_matlab_3D_array(EyR, outdims[2]); - free_cast_matlab_3D_array(EyI, outdims[2]); - free_cast_matlab_3D_array(EzR, outdims[2]); - free_cast_matlab_3D_array(EzI, outdims[2]); -} - -/** - * @brief Interpolate the magnetic field to the origin of the Yee cell - * - * @param[in] Hx_yee Steady state x component of magnetic field calculated at points in the Yee cell - * @param[in] Hy_yee Steady state y component of magnetic field calculated at points in the Yee cell - * @param[in] Hz_yee Steady state z component of magnetic field calculated at points in the Yee cell - * @param[out] Hx Steady state x component of magnetic field interpolated to Yee cell origin - * @param[out] Hy Steady state y component of magnetic field interpolated to Yee cell origin - * @param[out] Hz Steady state z component of magnetic field interpolated to Yee cell origin - * @param[in] I Number of elements in the i direction of the FDTD grid - * @param[in] J Number of elements in the j direction of the FDTD grid - * @param[in] K Number of elements in the k direction of the FDTD grid - * @param[in] i_l Least i index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest j index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest k index into the FDTD grid to evaluate the field at. Should be <= K-2 - */ -void interpolateFieldCentralH( double ***Hx_yee, double ***Hy_yee, double ***Hz_yee, - double ***Hx , double ***Hy , double ***Hz , - int I , int J , int K , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - int i,j,k; - double res1, res2, res3, res4; - - // throw runtime_error("Entering interpolateFieldCentralH\n"); - //first check that the limits of field extraction of within range - checkInterpolationPoints(i_l, i_u, j_l, j_u, k_l, k_u, I, J, K); - - if(j_u= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest j index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest k index into the FDTD grid to evaluate the field at. Should be <= K-2 - */ -void interpolateFieldCentralH_TE( double ***Hx_yee, double ***Hy_yee, double ***Hz_yee, - double ***Hx , double ***Hy , double ***Hz , - int I , int J , int K , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - int i,j,k; - double res1, res2, res3, res4; - - //first check that the limits of field extraction of within range - checkInterpolationPoints(i_l, i_u, j_l, j_u, k_l, k_u, I, J, K); - - for(k=k_l;k<=k_u;k++) - for(j=j_l;j<=j_u;j++) - for(i=i_l;i<=i_u;i++){ - /* - res1 = interp1(Hx_yee[k+1][j][i+1],Hx_yee[k+1][j][i],Hx_yee[k+1][j][i-1],Hx_yee[k+1][j][i-2]); - res2 = interp1(Hx_yee[k ][j][i+1],Hx_yee[k ][j][i],Hx_yee[k ][j][i-1],Hx_yee[k ][j][i-2]); - res3 = interp1(Hx_yee[k-1][j][i+1],Hx_yee[k-1][j][i],Hx_yee[k-1][j][i-1],Hx_yee[k-1][j][i-2]); - res4 = interp1(Hx_yee[k-2][j][i+1],Hx_yee[k-2][j][i],Hx_yee[k-2][j][i-1],Hx_yee[k-2][j][i-2]); - Hx[k-k_l][j-j_l][i-i_l] = interp1(res1, res2, res3, res4); - - res1 = interp1(Hy_yee[k+1][j+1][i],Hy_yee[k+1][j][i],Hy_yee[k+1][j-1][i],Hy_yee[k+1][j-2][i]); - res2 = interp1(Hy_yee[k ][j+1][i],Hy_yee[k ][j][i],Hy_yee[k ][j-1][i],Hy_yee[k ][j-2][i]); - res3 = interp1(Hy_yee[k-1][j+1][i],Hy_yee[k-1][j][i],Hy_yee[k-1][j-1][i],Hy_yee[k-1][j-2][i]); - res4 = interp1(Hy_yee[k-2][j+1][i],Hy_yee[k-2][j][i],Hy_yee[k-2][j-1][i],Hy_yee[k-2][j-2][i]); - Hy[k-k_l][j-j_l][i-i_l] = interp1(res1, res2, res3, res4); - - res1 = interp1(Hz_yee[k][j+1][i+1],Hz_yee[k][j+1][i],Hz_yee[k][j+1][i-1],Hz_yee[k][j+1][i-2]); - res2 = interp1(Hz_yee[k][j ][i+1],Hz_yee[k][j ][i],Hz_yee[k][j ][i-1],Hz_yee[k][j ][i-2]); - res3 = interp1(Hz_yee[k][j-1][i+1],Hz_yee[k][j-1][i],Hz_yee[k][j-1][i-1],Hz_yee[k][j-1][i-2]); - res4 = interp1(Hz_yee[k][j-2][i+1],Hz_yee[k][j-2][i],Hz_yee[k][j-2][i-1],Hz_yee[k][j-2][i-2]); - Hz[k-k_l][j-j_l][i-i_l] = interp1(res1, res2, res3, res4); - */ - - Hx[k-k_l][j-j_l][i-i_l] = 0.; - - Hy[k-k_l][j-j_l][i-i_l] = 0.; - - res1 = interp1(Hz_yee[k][j+1][i+1],Hz_yee[k][j+1][i],Hz_yee[k][j+1][i-1],Hz_yee[k][j+1][i-2]); - res2 = interp1(Hz_yee[k][j ][i+1],Hz_yee[k][j ][i],Hz_yee[k][j ][i-1],Hz_yee[k][j ][i-2]); - res3 = interp1(Hz_yee[k][j-1][i+1],Hz_yee[k][j-1][i],Hz_yee[k][j-1][i-1],Hz_yee[k][j-1][i-2]); - res4 = interp1(Hz_yee[k][j-2][i+1],Hz_yee[k][j-2][i],Hz_yee[k][j-2][i-1],Hz_yee[k][j-2][i-2]); - Hz[k-k_l][j-j_l][i-i_l] = interp1(res1, res2, res3, res4); - } -} - -/** - * @brief Interpolate the TE magnetic field (x,y components) to the origin of the Yee cell - * - * @param[in] Hx_yee Steady state x component of magnetic field calculated at points in the Yee cell - * @param[in] Hy_yee Steady state y component of magnetic field calculated at points in the Yee cell - * @param[in] Hz_yee Steady state z component of magnetic field calculated at points in the Yee cell - * @param[out] Hx Steady state x component of magnetic field interpolated to Yee cell origin - * @param[out] Hy Steady state y component of magnetic field interpolated to Yee cell origin - * @param[out] Hz Steady state z component of magnetic field interpolated to Yee cell origin - * @param[in] I Number of elements in the i direction of the FDTD grid - * @param[in] J Number of elements in the j direction of the FDTD grid - * @param[in] K Number of elements in the k direction of the FDTD grid - * @param[in] i_l Least i index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest j index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest k index into the FDTD grid to evaluate the field at. Should be <= K-2 - */ -void interpolateFieldCentralH_TM( double ***Hx_yee, double ***Hy_yee, double ***Hz_yee, - double ***Hx , double ***Hy , double ***Hz , - int I , int J , int K , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - int i,j,k; - - //first check that the limits of field extraction of within range - checkInterpolationPoints(i_l, i_u, j_l, j_u, k_l, k_u, I, J, K); - - for(k=k_l;k<=k_u;k++) - for(j=j_l;j<=j_u;j++) - for(i=i_l;i<=i_u;i++){ - /* - res1 = interp1(Hx_yee[k+1][j][i+1],Hx_yee[k+1][j][i],Hx_yee[k+1][j][i-1],Hx_yee[k+1][j][i-2]); - res2 = interp1(Hx_yee[k ][j][i+1],Hx_yee[k ][j][i],Hx_yee[k ][j][i-1],Hx_yee[k ][j][i-2]); - res3 = interp1(Hx_yee[k-1][j][i+1],Hx_yee[k-1][j][i],Hx_yee[k-1][j][i-1],Hx_yee[k-1][j][i-2]); - res4 = interp1(Hx_yee[k-2][j][i+1],Hx_yee[k-2][j][i],Hx_yee[k-2][j][i-1],Hx_yee[k-2][j][i-2]); - Hx[k-k_l][j-j_l][i-i_l] = interp1(res1, res2, res3, res4); - - res1 = interp1(Hy_yee[k+1][j+1][i],Hy_yee[k+1][j][i],Hy_yee[k+1][j-1][i],Hy_yee[k+1][j-2][i]); - res2 = interp1(Hy_yee[k ][j+1][i],Hy_yee[k ][j][i],Hy_yee[k ][j-1][i],Hy_yee[k ][j-2][i]); - res3 = interp1(Hy_yee[k-1][j+1][i],Hy_yee[k-1][j][i],Hy_yee[k-1][j-1][i],Hy_yee[k-1][j-2][i]); - res4 = interp1(Hy_yee[k-2][j+1][i],Hy_yee[k-2][j][i],Hy_yee[k-2][j-1][i],Hy_yee[k-2][j-2][i]); - Hy[k-k_l][j-j_l][i-i_l] = interp1(res1, res2, res3, res4); - - res1 = interp1(Hz_yee[k][j+1][i+1],Hz_yee[k][j+1][i],Hz_yee[k][j+1][i-1],Hz_yee[k][j+1][i-2]); - res2 = interp1(Hz_yee[k][j ][i+1],Hz_yee[k][j ][i],Hz_yee[k][j ][i-1],Hz_yee[k][j ][i-2]); - res3 = interp1(Hz_yee[k][j-1][i+1],Hz_yee[k][j-1][i],Hz_yee[k][j-1][i-1],Hz_yee[k][j-1][i-2]); - res4 = interp1(Hz_yee[k][j-2][i+1],Hz_yee[k][j-2][i],Hz_yee[k][j-2][i-1],Hz_yee[k][j-2][i-2]); - Hz[k-k_l][j-j_l][i-i_l] = interp1(res1, res2, res3, res4); - */ - - Hx[k-k_l][j-j_l][i-i_l] = interp1(Hx_yee[k][j-2][i], Hx_yee[k][j-1][i], Hx_yee[k][j][i], Hx_yee[k][j+1][i]); - - Hy[k-k_l][j-j_l][i-i_l] = interp1(Hy_yee[k][j][i-2], Hy_yee[k][j][i-1], Hy_yee[k][j][i], Hy_yee[k][j][i+1]); - - Hz[k-k_l][j-j_l][i-i_l] = 0.; - } -} - -/** - * @brief Interpolate the magnetic field to the origin of the Yee cell - MATLAB interface - * - * @param[in] Hx_yee Input Yee cell Hx field matrix - * @param[in] Hy_yee Input Yee cell Hy field matrix - * @param[in] Hz_yee Input Yee cell Hz field matrix - * @param[out] Hx Output Yee cell Hx field matrix - * @param[out] Hy Output Yee cell Hy field matrix - * @param[out] Hz Output Yee cell Hz field matrix - * @param[in] i_l Least i index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= K-2 - * - * I, J, K are the number of elements in the i, j, k directions of the FDTD grid, respectively. - */ -void mxInterpolateFieldCentralH( mxArray *Hx_yee , mxArray *Hy_yee , mxArray *Hz_yee, - mxArray **Hx , mxArray **Hy , mxArray **Hz , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - double ***HxR, ***HxI, ***HyR, ***HyI, ***HzR, ***HzI,***Hx_yee_R, ***Hx_yee_I, ***Hy_yee_R, ***Hy_yee_I, ***Hz_yee_R, ***Hz_yee_I; - const int *indims; - int outdims[3], ndims; - - if( mxGetNumberOfDimensions(Hx_yee) < 3){ - throw runtime_error("Error in mxInterpolateFieldCentralH, Ex_yee does not have 3 dimensions\n"); - } - - indims = (int *)mxGetDimensions(Hx_yee); - //assume that all matrices have the same dimensions - if( !mxIsComplex(Hx_yee) ){ - mexErrMsgTxt("Ex_yee is not complex"); - } - if( !mxIsComplex(Hy_yee) ){ - mexErrMsgTxt("Ey_yee is not complex"); - } - if( !mxIsComplex(Hz_yee) ){ - mexErrMsgTxt("Ez_yee is not complex"); - } - - - Hx_yee_R = cast_matlab_3D_array(mxGetPr(Hx_yee), indims[0], indims[1], indims[2]); - Hx_yee_I = cast_matlab_3D_array(mxGetPi(Hx_yee), indims[0], indims[1], indims[2]); - - Hy_yee_R = cast_matlab_3D_array(mxGetPr(Hy_yee), indims[0], indims[1], indims[2]); - Hy_yee_I = cast_matlab_3D_array(mxGetPi(Hy_yee), indims[0], indims[1], indims[2]); - - Hz_yee_R = cast_matlab_3D_array(mxGetPr(Hz_yee), indims[0], indims[1], indims[2]); - Hz_yee_I = cast_matlab_3D_array(mxGetPi(Hz_yee), indims[0], indims[1], indims[2]); - - //now construct the output matrices - ndims = 3; - - outdims[0] = i_u - i_l + 1; - outdims[1] = j_u - j_l + 1; - if(outdims[1]<1) - outdims[1]=1; - outdims[2] = k_u - k_l + 1; - - *Hx = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Hy = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Hz = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - - HxR = cast_matlab_3D_array(mxGetPr(*Hx), outdims[0], outdims[1], outdims[2]); - HxI = cast_matlab_3D_array(mxGetPi(*Hx), outdims[0], outdims[1], outdims[2]); - - HyR = cast_matlab_3D_array(mxGetPr(*Hy), outdims[0], outdims[1], outdims[2]); - HyI = cast_matlab_3D_array(mxGetPi(*Hy), outdims[0], outdims[1], outdims[2]); - - HzR = cast_matlab_3D_array(mxGetPr(*Hz), outdims[0], outdims[1], outdims[2]); - HzI = cast_matlab_3D_array(mxGetPi(*Hz), outdims[0], outdims[1], outdims[2]); - - //now finally ready for interpolation - interpolateFieldCentralH( Hx_yee_I, Hy_yee_I, Hz_yee_I, - HxI , HyI , HzI , - indims[0] , indims[1] , indims[2] , - i_l, i_u, j_l, j_u, k_l, k_u); - - - interpolateFieldCentralH( Hx_yee_R, Hy_yee_R, Hz_yee_R, - HxR , HyR , HzR , - indims[0] , indims[1] , indims[2] , - i_l, i_u, j_l, j_u, k_l, k_u); - - //free the extra memory used by casting array - free_cast_matlab_3D_array(Hx_yee_R, indims[2]); - free_cast_matlab_3D_array(Hx_yee_I, indims[2]); - free_cast_matlab_3D_array(Hy_yee_R, indims[2]); - free_cast_matlab_3D_array(Hy_yee_I, indims[2]); - free_cast_matlab_3D_array(Hz_yee_R, indims[2]); - free_cast_matlab_3D_array(Hz_yee_I, indims[2]); - - free_cast_matlab_3D_array(HxR, outdims[2]); - free_cast_matlab_3D_array(HxI, outdims[2]); - free_cast_matlab_3D_array(HyR, outdims[2]); - free_cast_matlab_3D_array(HyI, outdims[2]); - free_cast_matlab_3D_array(HzR, outdims[2]); - free_cast_matlab_3D_array(HzI, outdims[2]); -} - -/** - * @brief Interpolate the TE magnetic field (z component) to the origin of the Yee cell - MATLAB interface - * - * @param[in] Hx_yee Input Yee cell Hx field matrix - * @param[in] Hy_yee Input Yee cell Hy field matrix - * @param[in] Hz_yee Input Yee cell Hz field matrix - * @param[out] Hx Output Yee cell Hx field matrix - * @param[out] Hy Output Yee cell Hy field matrix - * @param[out] Hz Output Yee cell Hz field matrix - * @param[in] i_l Least i index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= K-2 - * - * I, J, K are the number of elements in the i, j, k directions of the FDTD grid, respectively. - */ -void mxInterpolateFieldCentralH_TE( mxArray *Hx_yee , mxArray *Hy_yee , mxArray *Hz_yee, - mxArray **Hx , mxArray **Hy , mxArray **Hz , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - double ***HxR, ***HxI, ***HyR, ***HyI, ***HzR, ***HzI,***Hx_yee_R, ***Hx_yee_I, ***Hy_yee_R, ***Hy_yee_I, ***Hz_yee_R, ***Hz_yee_I; - const int *indims; - int outdims[3], ndims; - - if( mxGetNumberOfDimensions(Hx_yee) != 2){ - throw runtime_error("Error in mxInterpolateFieldCentralH_TE, Ex_yee does not have 3 dimensions\n"); - - - } - - indims = (int *)mxGetDimensions(Hx_yee); - //assume that all matrices have the same dimensions - if( !mxIsComplex(Hx_yee) ){ - mexErrMsgTxt("Ex_yee is not complex"); - } - if( !mxIsComplex(Hy_yee) ){ - mexErrMsgTxt("Ey_yee is not complex"); - } - if( !mxIsComplex(Hz_yee) ){ - mexErrMsgTxt("Ez_yee is not complex"); - } - - - Hx_yee_R = cast_matlab_3D_array(mxGetPr(Hx_yee), indims[0], indims[1], 0); - Hx_yee_I = cast_matlab_3D_array(mxGetPi(Hx_yee), indims[0], indims[1], 0); - - Hy_yee_R = cast_matlab_3D_array(mxGetPr(Hy_yee), indims[0], indims[1], 0); - Hy_yee_I = cast_matlab_3D_array(mxGetPi(Hy_yee), indims[0], indims[1], 0); - - Hz_yee_R = cast_matlab_3D_array(mxGetPr(Hz_yee), indims[0], indims[1], 0); - Hz_yee_I = cast_matlab_3D_array(mxGetPi(Hz_yee), indims[0], indims[1], 0); - - //now construct the output matrices - ndims = 3; - - outdims[0] = i_u - i_l + 1; - outdims[1] = j_u - j_l + 1; - outdims[2] = 1; - - *Hx = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Hy = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Hz = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - - HxR = cast_matlab_3D_array(mxGetPr(*Hx), outdims[0], outdims[1], outdims[2]); - HxI = cast_matlab_3D_array(mxGetPi(*Hx), outdims[0], outdims[1], outdims[2]); - - HyR = cast_matlab_3D_array(mxGetPr(*Hy), outdims[0], outdims[1], outdims[2]); - HyI = cast_matlab_3D_array(mxGetPi(*Hy), outdims[0], outdims[1], outdims[2]); - - HzR = cast_matlab_3D_array(mxGetPr(*Hz), outdims[0], outdims[1], outdims[2]); - HzI = cast_matlab_3D_array(mxGetPi(*Hz), outdims[0], outdims[1], outdims[2]); - - //now finally ready for interpolation - interpolateFieldCentralH_TE( Hx_yee_I, Hy_yee_I, Hz_yee_I, - HxI , HyI , HzI , - indims[0] , indims[1] , indims[2] , - i_l, i_u, j_l, j_u, k_l, k_u); - - - interpolateFieldCentralH_TE( Hx_yee_R, Hy_yee_R, Hz_yee_R, - HxR , HyR , HzR , - indims[0] , indims[1] , indims[2] , - i_l, i_u, j_l, j_u, k_l, k_u); - - //free the extra memory used by casting array - free_cast_matlab_3D_array(Hx_yee_R, 0); - free_cast_matlab_3D_array(Hx_yee_I, 0); - free_cast_matlab_3D_array(Hy_yee_R, 0); - free_cast_matlab_3D_array(Hy_yee_I, 0); - free_cast_matlab_3D_array(Hz_yee_R, 0); - free_cast_matlab_3D_array(Hz_yee_I, 0); - - free_cast_matlab_3D_array(HxR, outdims[2]); - free_cast_matlab_3D_array(HxI, outdims[2]); - free_cast_matlab_3D_array(HyR, outdims[2]); - free_cast_matlab_3D_array(HyI, outdims[2]); - free_cast_matlab_3D_array(HzR, outdims[2]); - free_cast_matlab_3D_array(HzI, outdims[2]); -} - -/** - * @brief Interpolate the TM magnetic field (x,y components) to the origin of the Yee cell - MATLAB interface - * - * @param[in] Hx_yee Input Yee cell Hx field matrix - * @param[in] Hy_yee Input Yee cell Hy field matrix - * @param[in] Hz_yee Input Yee cell Hz field matrix - * @param[out] Hx Output Yee cell Hx field matrix - * @param[out] Hy Output Yee cell Hy field matrix - * @param[out] Hz Output Yee cell Hz field matrix - * @param[in] i_l Least i index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] i_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= I-2 - * @param[in] j_l Least j index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] j_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= J-2 - * @param[in] k_l Least k index into the FDTD grid to evaluate the field at. Should be >= 2 - * @param[in] k_u Greatest i index into the FDTD grid to evaluate the field at. Should be <= K-2 - * - * I, J, K are the number of elements in the i, j, k directions of the FDTD grid, respectively. - */ -void mxInterpolateFieldCentralH_TM( mxArray *Hx_yee , mxArray *Hy_yee , mxArray *Hz_yee, - mxArray **Hx , mxArray **Hy , mxArray **Hz , - int i_l, int i_u, int j_l, int j_u, int k_l, int k_u){ - - double ***HxR, ***HxI, ***HyR, ***HyI, ***HzR, ***HzI,***Hx_yee_R, ***Hx_yee_I, ***Hy_yee_R, ***Hy_yee_I, ***Hz_yee_R, ***Hz_yee_I; - const int *indims; - int outdims[3], ndims; - - if( mxGetNumberOfDimensions(Hx_yee) != 2){ - throw runtime_error("Error in mxInterpolateFieldCentralH_TM, Ex_yee does not have 3 dimensions\n"); - } - - indims = (int *)mxGetDimensions(Hx_yee); - //assume that all matrices have the same dimensions - if( !mxIsComplex(Hx_yee) ){ - mexErrMsgTxt("Ex_yee is not complex"); - } - if( !mxIsComplex(Hy_yee) ){ - mexErrMsgTxt("Ey_yee is not complex"); - } - if( !mxIsComplex(Hz_yee) ){ - mexErrMsgTxt("Ez_yee is not complex"); - } - - - Hx_yee_R = cast_matlab_3D_array(mxGetPr(Hx_yee), indims[0], indims[1], 0); - Hx_yee_I = cast_matlab_3D_array(mxGetPi(Hx_yee), indims[0], indims[1], 0); - - Hy_yee_R = cast_matlab_3D_array(mxGetPr(Hy_yee), indims[0], indims[1], 0); - Hy_yee_I = cast_matlab_3D_array(mxGetPi(Hy_yee), indims[0], indims[1], 0); - - Hz_yee_R = cast_matlab_3D_array(mxGetPr(Hz_yee), indims[0], indims[1], 0); - Hz_yee_I = cast_matlab_3D_array(mxGetPi(Hz_yee), indims[0], indims[1], 0); - - //now construct the output matrices - ndims = 3; - - outdims[0] = i_u - i_l + 1; - outdims[1] = j_u - j_l + 1; - outdims[2] = 1; - - *Hx = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Hy = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - *Hz = mxCreateNumericArray( ndims, (const mwSize *)outdims, mxDOUBLE_CLASS, mxCOMPLEX); - - HxR = cast_matlab_3D_array(mxGetPr(*Hx), outdims[0], outdims[1], outdims[2]); - HxI = cast_matlab_3D_array(mxGetPi(*Hx), outdims[0], outdims[1], outdims[2]); - - HyR = cast_matlab_3D_array(mxGetPr(*Hy), outdims[0], outdims[1], outdims[2]); - HyI = cast_matlab_3D_array(mxGetPi(*Hy), outdims[0], outdims[1], outdims[2]); - - HzR = cast_matlab_3D_array(mxGetPr(*Hz), outdims[0], outdims[1], outdims[2]); - HzI = cast_matlab_3D_array(mxGetPi(*Hz), outdims[0], outdims[1], outdims[2]); - - //now finally ready for interpolation - interpolateFieldCentralH_TM( Hx_yee_I, Hy_yee_I, Hz_yee_I, - HxI , HyI , HzI , - indims[0] , indims[1] , indims[2] , - i_l, i_u, j_l, j_u, k_l, k_u); - - - interpolateFieldCentralH_TM( Hx_yee_R, Hy_yee_R, Hz_yee_R, - HxR , HyR , HzR , - indims[0] , indims[1] , indims[2] , - i_l, i_u, j_l, j_u, k_l, k_u); - - //free the extra memory used by casting array - free_cast_matlab_3D_array(Hx_yee_R, 0); - free_cast_matlab_3D_array(Hx_yee_I, 0); - free_cast_matlab_3D_array(Hy_yee_R, 0); - free_cast_matlab_3D_array(Hy_yee_I, 0); - free_cast_matlab_3D_array(Hz_yee_R, 0); - free_cast_matlab_3D_array(Hz_yee_I, 0); - - free_cast_matlab_3D_array(HxR, outdims[2]); - free_cast_matlab_3D_array(HxI, outdims[2]); - free_cast_matlab_3D_array(HyR, outdims[2]); - free_cast_matlab_3D_array(HyI, outdims[2]); - free_cast_matlab_3D_array(HzR, outdims[2]); - free_cast_matlab_3D_array(HzI, outdims[2]); -} - -/*Interpolate the electric field to the origin of the Yee cell in the time domain - * - *[in]Hxy to Hzy are the split components of the Yee cell - *[in](i,j,k) denotes the particular Yee cell we are interested in - * - *[out]Hx to Hz are the field components interpolated to the origin of the Yee cell - */ -void interpolateTimeDomainFieldCentralH( SplitFieldComponent& Hxy, SplitFieldComponent& Hxz, SplitFieldComponent& Hyx, SplitFieldComponent& Hyz, SplitFieldComponent& Hzx, SplitFieldComponent& Hzy, - int i, int j, int k, - double *Hx, double *Hy, double *Hz){ - - - double res1, res2, res3, res4; - - res1 = interp1(Hxy[k+1][j+1][i]+Hxz[k+1][j+1][i],Hxy[k+1][j][i]+Hxz[k+1][j][i],Hxy[k+1][j-1][i]+Hxz[k+1][j-1][i],Hxy[k+1][j-2][i]+Hxz[k+1][j-2][i]); - res2 = interp1(Hxy[k ][j+1][i]+Hxz[k ][j+1][i],Hxy[k ][j][i]+Hxz[k ][j][i],Hxy[k ][j-1][i]+Hxz[k ][j-1][i],Hxy[k ][j-2][i]+Hxz[k ][j-2][i]); - res3 = interp1(Hxy[k-1][j+1][i]+Hxz[k-1][j+1][i],Hxy[k-1][j][i]+Hxz[k-1][j][i],Hxy[k-1][j-1][i]+Hxz[k-1][j-1][i],Hxy[k-1][j-2][i]+Hxz[k-1][j-2][i]); - res4 = interp1(Hxy[k-2][j+1][i]+Hxz[k-2][j+1][i],Hxy[k-2][j][i]+Hxz[k-2][j][i],Hxy[k-2][j-1][i]+Hxz[k-2][j-1][i],Hxy[k-2][j-2][i]+Hxz[k-2][j-2][i]); - *Hx = interp1(res1, res2, res3, res4); - - res1 = interp1(Hyx[k+1][j][i+1]+Hyz[k+1][j][i+1],Hyx[k+1][j][i]+Hyz[k+1][j][i],Hyx[k+1][j][i-1]+Hyz[k+1][j][i-1],Hyx[k+1][j][i-2]+Hyz[k+1][j][i-2]); - res2 = interp1(Hyx[k ][j][i+1]+Hyz[k ][j][i+1],Hyx[k ][j][i]+Hyz[k ][j][i],Hyx[k ][j][i-1]+Hyz[k ][j][i-1],Hyx[k ][j][i-2]+Hyz[k ][j][i-2]); - res3 = interp1(Hyx[k-1][j][i+1]+Hyz[k-1][j][i+1],Hyx[k-1][j][i]+Hyz[k-1][j][i],Hyx[k-1][j][i-1]+Hyz[k-1][j][i-1],Hyx[k-1][j][i-2]+Hyz[k-1][j][i-2]); - res4 = interp1(Hyx[k-2][j][i+1]+Hyz[k-2][j][i+1],Hyx[k-2][j][i]+Hyz[k-2][j][i],Hyx[k-2][j][i-1]+Hyz[k-2][j][i-1],Hyx[k-2][j][i-2]+Hyz[k-2][j][i-2]); - *Hy = interp1(res1, res2, res3, res4); - - res1 = interp1(Hzx[k][j+1][i+1]+Hzy[k][j+1][i+1],Hzx[k][j+1][i]+Hzy[k][j+1][i],Hzx[k][j+1][i-1]+Hzy[k][j+1][i-1],Hzx[k][j+1][i-2]+Hzy[k][j+1][i-2]); - res2 = interp1(Hzx[k][j ][i+1]+Hzy[k][j ][i+1],Hzx[k][j ][i]+Hzy[k][j ][i],Hzx[k][j ][i-1]+Hzy[k][j ][i-1],Hzx[k][j ][i-2]+Hzy[k][j ][i-2]); - res3 = interp1(Hzx[k][j-1][i+1]+Hzy[k][j-1][i+1],Hzx[k][j-1][i]+Hzy[k][j-1][i],Hzx[k][j-1][i-1]+Hzy[k][j-1][i-1],Hzx[k][j-1][i-2]+Hzy[k][j-1][i-2]); - res4 = interp1(Hzx[k][j-2][i+1]+Hzy[k][j-2][i+1],Hzx[k][j-2][i]+Hzy[k][j-2][i],Hzx[k][j-2][i-1]+Hzy[k][j-2][i-1],Hzx[k][j-2][i-2]+Hzy[k][j-2][i-2]); - *Hz = interp1(res1, res2, res3, res4); -} - -/*Interpolate the electric field to the origin of the Yee cell in the time domain using band limited interpolation - * - *[in]Hxy to Hzy are the split components of the Yee cell - *[in](i,j,k) denotes the particular Yee cell we are interested in - * - *[out]Hx to Hz are the field components interpolated to the origin of the Yee cell - */ -void interpolateTimeDomainFieldCentralHBandLimited( SplitFieldComponent& Hxy, SplitFieldComponent& Hxz, SplitFieldComponent& Hyx, SplitFieldComponent& Hyz, SplitFieldComponent& Hzx, SplitFieldComponent& Hzy, - int i, int j, int k, - double *Hx, double *Hy, double *Hz){ - - /*Array for performing bandwidth limited interpolation obtained using Matlab's interp function*/ - const int Nbvec = 8; - const double bvec[Nbvec] = {-0.006777513830606,0.039457774230186,-0.142658093428622,0.609836360661632,0.609836360661632,-0.142658093428622,0.039457774230186,-0.006777513830606}; - *Hx = 0.; - *Hy = 0.; - *Hz = 0.; - - double hx = 0; - double hy = 0; - double hz = 0; - - - for(int ind1=0;ind1 nI - 2) - { - throw runtime_error("Interpolation error: i_u too large"); - } - else if (j_l < 2) - { - throw runtime_error("Interpolation error: j_l too small"); - } - else if (j_u > nJ - 2) - { - throw runtime_error("Interpolation error: j_u too large"); - } - else if (k_l < 2) - { - throw runtime_error("Interpolation error: k_l too small"); - } - else if (k_u > nK - 2) - { - throw runtime_error("Interpolation error: k_u too large"); - } -} - -// THESE WILL BE THE RELEVANT FUNCTIONS ONCE INTERPOLATION SCHEMES ARE COMPLETE AND TESTED - InterpolationScheme::InterpolationScheme(scheme_value val) { // set the value field @@ -241,64 +185,54 @@ bool InterpolationScheme::is_better_than(const InterpolationScheme s) const { return (priority > s.get_priority()); } -const InterpolationScheme &best_scheme(int cells_in_direction, int cell_id) { - - // interpolation is impossible with fewer than 4 datapoints in a dimension - // 4 datapoints => 3 cells in a given direction - if (cells_in_direction < 3) { - throw out_of_range("Error: domain axis has < 3 cells (< 4 datapoints), cubic and bandlimited interpolation impossible.\n"); - } - // Yee cell with index <0 doesn't allow for interpolation (this is the PML) - else if (cell_id < 0) { - throw out_of_range("Error: Yee cell index < 0 requested.\n"); - } - // Yee cell with index >= cells_in_direction doesn't exist - // Again, we can in theory determine the value "here" using BL7, however cubic interpolation cannot do this, and so current code functionality is to throw an error here. - else if (cell_id >= cells_in_direction) { - throw out_of_range("Error: Yee cell index beyond maximum number of Yee cells requested.\n"); - } - else if (cells_in_direction < 7) { - // we do not have enough cells to use bandlimited interpolation, but can use cubic - // cell_id = 0 requires us to use CBFst, - // cell_id = 1,2,...,cells_in_direction-2 requires us to use CBMid, - // cell_id = cells_in_direction-1 requires us to use CBLast - if (cell_id == 0) { - return CBFst; - } - else if (cell_id == cells_in_direction-1) { - return CBLst; - } - else { - return CBMid; - } +const InterpolationScheme &best_scheme(int datapts_in_direction, int interpolation_position, + PreferredInterpolationMethods interpolation_methods) { + // interpolation is impossible with fewer than 4 datapoints in a dimension + if (datapts_in_direction < 4) { + throw out_of_range("Error: domain axis has <4 datapoints, cubic and bandlimited interpolation " + "impossible.\n"); + } else if ((datapts_in_direction < 8) || (interpolation_methods == PreferredInterpolationMethods::Cubic)) { + // we are restricted to cubic interpolation + if (interpolation_position <= 0 || interpolation_position >= datapts_in_direction) { + throw out_of_range("Error: Cubic interpolation impossible to position " + + to_string(interpolation_position) + "\n"); + } else { + // determine cubic interpolation scheme to use + if (interpolation_position == 1) { + return CBFst; + } else if (interpolation_position == datapts_in_direction - 1) { + return CBLst; + } else { + return CBMid; + } } - else { - // we can apply bandlimited interpolation. - // if there are 4 datapoints either side, we want to use BL3 - if ((cell_id > 2) && (cell_id < cells_in_direction - 3)) { - return BL3; - } - else if (cell_id == 0) { - return BL0; - } - else if (cell_id == 1) { - return BL1; - } - else if (cell_id == 2) { - return BL2; - } - else if (cell_id == cells_in_direction - 3) { - return BL4; - } - else if (cell_id == cells_in_direction - 2) { - return BL5; - } - else if (cell_id == cells_in_direction - 1) { - return BL6; - } - else { - // we somehow got to here, but we should have covered all possible bases above. Return an error - throw runtime_error("Error: could not identify scheme despite appropriate Yee cell index and number of cells.\n"); - } + } else if (interpolation_position < 0 || interpolation_position > datapts_in_direction) { + // cannot interpolate to here using BLi, throw error + throw out_of_range("Error: BLi interpolation impossible to position " + + to_string(interpolation_position) + "\n"); + } else { + // safe to use BLi, figure out which scheme we need + if (interpolation_position == 0) { + return BL_TO_CELL_0; + } else if (interpolation_position == datapts_in_direction) { + return BL7; + } else if (interpolation_position == 1) { + return BL0; + } else if (interpolation_position == 2) { + return BL1; + } else if (interpolation_position == 3) { + return BL2; + } else if (interpolation_position == datapts_in_direction - 3) { + return BL4; + } else if (interpolation_position == datapts_in_direction - 2) { + return BL5; + } else if (interpolation_position == datapts_in_direction - 1) { + return BL6; + } else { + // we have 4 datapoints either side of where we want to interpolate to + return BL3; } + } + // if we get to here we have, somehow, not returned a scheme. Raise an error + throw runtime_error("Error: could not identify scheme for unknown reasons, diagnose further.\n"); } diff --git a/tdms/src/iterator.cpp b/tdms/src/iterator.cpp index d69ff8376..00b5c6612 100644 --- a/tdms/src/iterator.cpp +++ b/tdms/src/iterator.cpp @@ -12,7 +12,6 @@ #include "array_init.h" #include "globals.h" #include "interface.h" -#include "interpolate.h" #include "matlabio.h" #include "shapes.h" #include "source.h" @@ -29,7 +28,7 @@ using namespace tdms_phys_constants; #define TIME_MAIN_LOOP true //threshold used to terminate the steady state iterations #define TOL 1e-6 -//parameter to control the with of the ramp when introducing the waveform in steady state mode +//parameter to control the PreferredInterpolationMethodswith of the ramp when introducing the waveform in steady state mode #define ramp_width 4. /*This mex function will take in the following arguments and perform the @@ -233,23 +232,35 @@ using namespace tdms_phys_constants; campssample.components - numerical array of up to six elements which defines which field components will be sampled, 1 means Ex, 2 Ey etc. */ -void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[], SolverMethod method) { +void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[], + SolverMethod solver_method, + PreferredInterpolationMethods preferred_interpolation_methods) { - if (method == SolverMethod::FiniteDifference) { + if (solver_method == SolverMethod::FiniteDifference) { spdlog::info("Using finite-difference method (FDTD)"); } else { spdlog::info("Using pseudospectral method (PSTD)"); } + if (preferred_interpolation_methods == PreferredInterpolationMethods::BandLimited) { + spdlog::info("Using band-limited interpolation where possible"); + } else { + spdlog::info("Restricting to cubic interpolation"); + } auto params = SimulationParameters(); auto E_s = ElectricSplitField(); + E_s.set_preferred_interpolation_methods(preferred_interpolation_methods); auto H_s = MagneticSplitField(); + H_s.set_preferred_interpolation_methods(preferred_interpolation_methods); auto J_s = CurrentDensitySplitField(); auto E = ElectricField(); + E.set_preferred_interpolation_methods(preferred_interpolation_methods); auto H = MagneticField(); + H.set_preferred_interpolation_methods(preferred_interpolation_methods); auto E_copy = ElectricField(); // Used to check convergence with E - E_copy + E_copy.set_preferred_interpolation_methods(preferred_interpolation_methods); // We never actually interpolate this field, but adding this just in case we later add functionality that depends on it double ***surface_EHr, ***surface_EHi; double rho; @@ -282,7 +293,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs int Ni_tdf = 0, Nk_tdf = 0; int skip_tdf = 1; - if (method == SolverMethod::FiniteDifference) skip_tdf = 6; + if (solver_method == SolverMethod::FiniteDifference) skip_tdf = 6; // PSTD storage (not used if FD) fftw_complex *dk_e_x, *dk_e_y, *dk_e_z, *dk_h_x, *dk_h_y, *dk_h_z; @@ -526,7 +537,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs CCoefficientMatrix ca_vec, cb_vec, cc_vec; EHVec eh_vec; - if (method == SolverMethod::PseudoSpectral) { + if (solver_method == SolverMethod::PseudoSpectral) { int max_IJK = E_s.max_IJK_tot(), n_threads = omp_get_max_threads(); ca_vec.allocate(n_threads, max_IJK + 1); cb_vec.allocate(n_threads, max_IJK + 1); @@ -558,7 +569,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs init_diff_shift_op(0.5, dk_h_x, N_h_x); init_diff_shift_op(0.5, dk_h_y, N_h_y); init_diff_shift_op(0.5, dk_h_z, N_h_z); - } // if (method == DerivativeMethod::PseudoSpectral) + } // if (solver_method == DerivativeMethod::PseudoSpectral) params.set_Np(f_ex_vec); @@ -1173,20 +1184,22 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs for (int kt = 0; kt < fieldsample.k.size(); kt++) for (int jt = 0; jt < fieldsample.j.size(); jt++) for (int it = 0; it < fieldsample.i.size(); it++) { - ////fprintf(stderr,"Pos fs 1\n"); - interpolateTimeDomainFieldCentralEBandLimited( - E_s.xy, E_s.xz, E_s.yx, E_s.yz, E_s.zx, E_s.zy, - fieldsample.i[it] + params.pml.Dxl - 1, - fieldsample.j[jt] + params.pml.Dyl - 1, - fieldsample.k[kt] + params.pml.Dzl - 1, &Ex_temp, &Ey_temp, &Ez_temp); - //fprintf(stderr,"Pos fs 2\n"); + CellCoordinate current_cell(fieldsample.i[it] + params.pml.Dxl - 1, + fieldsample.j[jt] + params.pml.Dyl - 1, + fieldsample.k[kt] + params.pml.Dzl - 1); + Ex_temp = E_s.interpolate_to_centre_of(AxialDirection::X, current_cell); + if (current_cell.j() != 0) { + Ey_temp = E_s.interpolate_to_centre_of(AxialDirection::Y, current_cell); + } else { + Ey_temp = E_s.yx[current_cell] + E_s.yz[current_cell]; + } + Ez_temp = E_s.interpolate_to_centre_of(AxialDirection::Z, current_cell); for (int nt = 0; nt < fieldsample.n.size(); nt++) fieldsample[nt][kt][jt][it] = fieldsample[nt][kt][jt][it] + pow(Ex_temp * Ex_temp + Ey_temp * Ey_temp + Ez_temp * Ez_temp, fieldsample.n[nt] / 2.) / params.Nt; - //fprintf(stderr,"%d %d %d %d -> %d %d %d (%d) %d [%d %d]\n",nt,kt,jt,it,(int)fieldsample_vecs.n[nt], (int)fieldsample_vecs.i[it] + params.pml.Dxl - 1, (int)fieldsample_vecs.j[jt] + params.pml.Dyl - 1, params.pml.Dyl,(int)fieldsample_vecs.k[kt] + params.pml.Dzl - 1 , Nsteps, (int)fieldsample_vecs.n[nt] - 2); } } } @@ -1356,8 +1369,8 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs Enp1 = 0.0; array_ind = 0; - if (params.dimension == THREE || params.dimension == TE) { - if (method == SolverMethod::FiniteDifference) { + if (params.dimension == THREE || params.dimension == Dimension::TRANSVERSE_ELECTRIC) { + if (solver_method == SolverMethod::FiniteDifference) { //FDTD, E_s.xy #pragma omp for for (k = 0; k < (K_tot + 1); k++) @@ -1624,7 +1637,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } } //PSTD, E_s.xy - } // if (method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) + } // if (solver_method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) /* if(is_disp){ i=36; @@ -1637,7 +1650,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs //fprintf(stderr,"Pos 04:\n"); //E_s.xz updates - if (method == SolverMethod::FiniteDifference) { + if (solver_method == SolverMethod::FiniteDifference) { #pragma omp for for (k = 1; k < K_tot; k++) for (j = 0; j < J_tot_p1_bound; j++) @@ -1908,11 +1921,11 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } } //PSTD, E_s.xz - } // if (method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) + } // if (solver_method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) //fprintf(stderr,"Pos 05:\n"); //E_s.yx updates - if (method == SolverMethod::FiniteDifference) { + if (solver_method == SolverMethod::FiniteDifference) { //FDTD, E_s.yx #pragma omp for for (k = 0; k < (K_tot + 1); k++) @@ -2166,11 +2179,11 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } } //PSTD, E_s.yx - }// if (method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) + }// if (solver_method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) //fprintf(stderr,"Pos 06:\n"); //E_s.yz updates - if (method == SolverMethod::FiniteDifference) { + if (solver_method == SolverMethod::FiniteDifference) { //FDTD, E_s.yz #pragma omp for for (k = 1; k < K_tot; k++) @@ -2416,12 +2429,12 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } } //PSTD, E_s.yz - }// if (method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) + }// if (solver_method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) } //if(params.dimension==THREE || params.dimension==TE) //fprintf(stderr,"Pos 07:\n"); - if (params.dimension == THREE || params.dimension == TE) { - if (method == SolverMethod::FiniteDifference) { + if (params.dimension == THREE || params.dimension == Dimension::TRANSVERSE_ELECTRIC) { + if (solver_method == SolverMethod::FiniteDifference) { #pragma omp for //E_s.zx updates for (k = 0; k < K_tot; k++) @@ -2679,7 +2692,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } } //PSTD, E_s.zx - }// if (method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) + }// if (solver_method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) }//(params.dimension==THREE || params.dimension==TE) else { #pragma omp for @@ -2765,8 +2778,8 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } } //fprintf(stderr,"Pos 08:\n"); - if (params.dimension == THREE || params.dimension == TE) { - if (method == SolverMethod::FiniteDifference) { + if (params.dimension == THREE || params.dimension == Dimension::TRANSVERSE_ELECTRIC) { + if (solver_method == SolverMethod::FiniteDifference) { //FDTD, E_s.zy #pragma omp for //E_s.zy updates @@ -3025,7 +3038,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } } //PSTD, E_s.zy - }// if (method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) + }// if (solver_method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) }//(params.dimension==THREE || params.dimension==TE) else { #pragma omp for @@ -3126,7 +3139,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs else array_ind = (I_tot + 1) * k + I0.index; - if (k < (K1.index) || params.dimension == TM) { + if (k < (K1.index) || params.dimension == Dimension::TRANSVERSE_MAGNETIC) { E_s.zx[k][j][I0.index] = E_s.zx[k][j][I0.index] - C.b.x[array_ind] * @@ -3173,7 +3186,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs else array_ind = (I_tot + 1) * k + I1.index; - if (k < (K1.index) || params.dimension == TM) { + if (k < (K1.index) || params.dimension == Dimension::TRANSVERSE_MAGNETIC) { E_s.zx[k][j][I1.index] = E_s.zx[k][j][I1.index] + C.b.x[array_ind] * @@ -3219,7 +3232,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs for (k = (K0.index); k <= (K1.index); k++) for (i = (I0.index); i <= (I1.index); i++) { if (J0.apply) {//Perform across J0 - if (k < (K1.index) || params.dimension == TM) { + if (k < (K1.index) || params.dimension == Dimension::TRANSVERSE_MAGNETIC) { if (!params.is_multilayer) array_ind = J0.index; else @@ -3271,7 +3284,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs else array_ind = (J_tot + 1) * k + J1.index; - if (k < (K1.index) || params.dimension == TM) { + if (k < (K1.index) || params.dimension == Dimension::TRANSVERSE_MAGNETIC) { E_s.zy[k][(J1.index)][i] = E_s.zy[k][(J1.index)][i] - C.b.y[array_ind] * @@ -3521,8 +3534,8 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs { n = omp_get_thread_num(); - if (params.dimension == THREE || params.dimension == TE) { - if (method == SolverMethod::FiniteDifference) { + if (params.dimension == THREE || params.dimension == Dimension::TRANSVERSE_ELECTRIC) { + if (solver_method == SolverMethod::FiniteDifference) { //FDTD, H_s.xz #pragma omp for //H_s.xz updates @@ -3597,9 +3610,9 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } //PSTD, H_s.xz - }// if (method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) + }// if (solver_method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) - if (method == SolverMethod::FiniteDifference) { + if (solver_method == SolverMethod::FiniteDifference) { //FDTD, H_s.xy #pragma omp for //H_s.xy updates @@ -3694,9 +3707,9 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs */ } //PSTD, H_s.xy - }// if (method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) + }// if (solver_method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) - if (method == SolverMethod::FiniteDifference) { + if (solver_method == SolverMethod::FiniteDifference) { //FDTD, H_s.yx #pragma omp for //H_s.yx updates @@ -3777,7 +3790,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs //PSTD, H_s.yx } - if (method == SolverMethod::FiniteDifference) { + if (solver_method == SolverMethod::FiniteDifference) { //FDTD, H_s.yz #pragma omp for //H_s.yz updates @@ -3868,7 +3881,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } } //PSTD, H_s.yz - }// if (method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) + }// if (solver_method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) } //(params.dimension==THREE || params.dimension==TE) else { @@ -3948,8 +3961,8 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } } - if (params.dimension == THREE || params.dimension == TE) { - if (method == SolverMethod::FiniteDifference) { + if (params.dimension == THREE || params.dimension == Dimension::TRANSVERSE_ELECTRIC) { + if (solver_method == SolverMethod::FiniteDifference) { //FDTD, H_s.zy #pragma omp for //H_s.zy update @@ -4027,10 +4040,10 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } } //PSTD, H_s.zy - }// if (method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) + }// if (solver_method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) - if (method == SolverMethod::FiniteDifference) { + if (solver_method == SolverMethod::FiniteDifference) { //FDTD, H_s.zx #pragma omp for //H_s.zx update @@ -4109,7 +4122,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs } } //PSTD, H_s.zx - }// if (method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) + }// if (solver_method == DerivativeMethod::FiniteDifference) (else PseudoSpectral) }//(params.dimension==THREE || params.dimension==TE) } //end parallel if (TIME_EXEC) { timer.click(); } @@ -4134,7 +4147,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs real(commonAmplitude * commonPhase * (Isource.real[k - (K0.index)][j - (J0.index)][0] + IMAGINARY_UNIT * Isource.imag[k - (K0.index)][j - (J0.index)][0])); - if (k < (K1.index) || params.dimension == TM) + if (k < (K1.index) || params.dimension == Dimension::TRANSVERSE_MAGNETIC) H_s.yx[k][j][(I0.index) - 1] = H_s.yx[k][j][(I0.index) - 1] - D.b.x[array_ind] * @@ -4155,7 +4168,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs real(commonAmplitude * commonPhase * (Isource.real[k - (K0.index)][j - (J0.index)][4] + IMAGINARY_UNIT * Isource.imag[k - (K0.index)][j - (J0.index)][4])); - if (k < (K1.index) || params.dimension == TM) + if (k < (K1.index) || params.dimension == Dimension::TRANSVERSE_MAGNETIC) H_s.yx[k][j][(I1.index)] = H_s.yx[k][j][(I1.index)] + D.b.x[array_ind] * @@ -4181,7 +4194,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs (Jsource.real[k - (K0.index)][i - (I0.index)][0] + IMAGINARY_UNIT * Jsource.imag[k - (K0.index)][i - (I0.index)][0])); - if (k < (K1.index) || params.dimension == TM) + if (k < (K1.index) || params.dimension == Dimension::TRANSVERSE_MAGNETIC) H_s.xy[k][(J0.index) - 1][i] = H_s.xy[k][(J0.index) - 1][i] + D.b.y[array_ind] * @@ -4202,7 +4215,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs real(commonAmplitude * commonPhase * (Jsource.real[k - (K0.index)][i - (I0.index)][4] + IMAGINARY_UNIT * Jsource.imag[k - (K0.index)][i - (I0.index)][4])); - if (k < (K1.index) || params.dimension == TM) + if (k < (K1.index) || params.dimension == Dimension::TRANSVERSE_MAGNETIC) H_s.xy[k][(J1.index)][i] = H_s.xy[k][(J1.index)][i] - D.b.y[array_ind] * @@ -4552,30 +4565,17 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs if (params.run_mode == RunMode::complete && params.exphasorsvolume) { //now interpolate over the extracted phasors if (params.dimension == THREE) { - fprintf(stderr, "mxInterpolateFieldCentralE: %d %d %d \n", E.I_tot - 2, - E.J_tot - 2, E.K_tot - 2); - //fprintf(stderr,"Pos 15_m1a\n"); - mxInterpolateFieldCentralE(plhs[0], plhs[1], plhs[2], &plhs[13], &plhs[14], &plhs[15], 2, - E.I_tot - 2, 2, E.J_tot - 2, 2, - E.K_tot - 2); - //fprintf(stderr,"Pos 15_m1b\n"); - - } else if (params.dimension == TE) - mxInterpolateFieldCentralE_TE(plhs[0], plhs[1], plhs[2], &plhs[13], &plhs[14], &plhs[15], 2, - E.I_tot - 2, 2, E.J_tot - 2, 0, 0); - else - mxInterpolateFieldCentralE_TM(plhs[0], plhs[1], plhs[2], &plhs[13], &plhs[14], &plhs[15], 2, - E.I_tot - 2, 2, E.J_tot - 2, 0, 0); - if (params.dimension == THREE) - mxInterpolateFieldCentralH(plhs[3], plhs[4], plhs[5], &plhs[16], &plhs[17], &plhs[18], 2, - E.I_tot - 2, 2, E.J_tot - 2, 2, - E.K_tot - 2); - else if (params.dimension == TE) - mxInterpolateFieldCentralH_TE(plhs[3], plhs[4], plhs[5], &plhs[16], &plhs[17], &plhs[18], 2, - E.I_tot - 2, 2, E.J_tot - 2, 0, 0); - else - mxInterpolateFieldCentralH_TM(plhs[3], plhs[4], plhs[5], &plhs[16], &plhs[17], &plhs[18], 2, - E.I_tot - 2, 2, E.J_tot - 2, 0, 0); + E.interpolate_over_range(&plhs[13], &plhs[14], &plhs[15], 2, E.I_tot - 2, 2, E.J_tot - 2, 2, + E.K_tot - 2, Dimension::THREE); + H.interpolate_over_range(&plhs[16], &plhs[17], &plhs[18], 2, H.I_tot - 2, 2, H.J_tot - 2, 2, + H.K_tot - 2, Dimension::THREE); + } else { + // either TE or TM, but interpolate_over_range will handle that for us. Only difference is the k_upper/lower values we pass... + E.interpolate_over_range(&plhs[13], &plhs[14], &plhs[15], 2, E.I_tot - 2, 2, E.J_tot - 2, 0, + 0, params.dimension); + H.interpolate_over_range(&plhs[16], &plhs[17], &plhs[18], 2, H.I_tot - 2, 2, H.J_tot - 2, 0, + 0, params.dimension); + } //fprintf(stderr,"Pos 15a\n"); //now set up the grid labels for the interpolated fields @@ -4737,7 +4737,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs free_cast_matlab_3D_array(materials, 0); /*Free the additional memory which was allocated to store integers which were passed as doubles*/ - if (method == SolverMethod::PseudoSpectral) { + if (solver_method == SolverMethod::PseudoSpectral) { fftw_free(dk_e_x); fftw_free(dk_e_y); fftw_free(dk_e_z); @@ -4873,50 +4873,37 @@ void extractPhasorsSurface(double **surface_EHr, double **surface_EHi, ElectricS { #pragma omp for for (vindex = 0; vindex < n_surface_vertices; vindex++) { - // fprintf(stderr,"vindex: %d: (%d %d %d)\n",vindex,surface_vertices[0][vindex],surface_vertices[1][vindex],surface_vertices[2][vindex]); - if (params.dimension == THREE) - if (J_tot == 0) { - interpolateTimeDomainFieldCentralE_2Dy( - E.xy, E.xz, E.yx, E.yz, E.zx, E.zy, surface_vertices[0][vindex], - surface_vertices[1][vindex], surface_vertices[2][vindex], &Ex, &Ey, &Ez); - } else if (params.interp_method == cubic) - interpolateTimeDomainFieldCentralE( - E.xy, E.xz, E.yx, E.yz, E.zx, E.zy, surface_vertices[0][vindex], - surface_vertices[1][vindex], surface_vertices[2][vindex], &Ex, &Ey, &Ez); - else - interpolateTimeDomainFieldCentralEBandLimited( - E.xy, E.xz, E.yx, E.yz, E.zx, E.zy, surface_vertices[0][vindex], - surface_vertices[1][vindex], surface_vertices[2][vindex], &Ex, &Ey, &Ez); - else if (params.dimension == TE) - interpolateTimeDomainFieldCentralE_TE( - E.xy, E.xz, E.yx, E.yz, E.zx, E.zy, surface_vertices[0][vindex], - surface_vertices[1][vindex], surface_vertices[2][vindex], &Ex, &Ey, &Ez); - else - interpolateTimeDomainFieldCentralE_TM( - E.xy, E.xz, E.yx, E.yz, E.zx, E.zy, surface_vertices[0][vindex], - surface_vertices[1][vindex], surface_vertices[2][vindex], &Ex, &Ey, &Ez); - // fprintf(stderr,"1st interp donezn"); - if (params.dimension == THREE) - if (J_tot == 0) { - interpolateTimeDomainFieldCentralH_2Dy( - H.xy, H.xz, H.yx, H.yz, H.zx, H.zy, surface_vertices[0][vindex], - surface_vertices[1][vindex], surface_vertices[2][vindex], &Hx, &Hy, &Hz); - } else if (params.interp_method == cubic) - interpolateTimeDomainFieldCentralH( - H.xy, H.xz, H.yx, H.yz, H.zx, H.zy, surface_vertices[0][vindex], - surface_vertices[1][vindex], surface_vertices[2][vindex], &Hx, &Hy, &Hz); - else - interpolateTimeDomainFieldCentralHBandLimited( - H.xy, H.xz, H.yx, H.yz, H.zx, H.zy, surface_vertices[0][vindex], - surface_vertices[1][vindex], surface_vertices[2][vindex], &Hx, &Hy, &Hz); - else if (params.dimension == TE) - interpolateTimeDomainFieldCentralH_TE( - H.xy, H.xz, H.yx, H.yz, H.zx, H.zy, surface_vertices[0][vindex], - surface_vertices[1][vindex], surface_vertices[2][vindex], &Hx, &Hy, &Hz); - else - interpolateTimeDomainFieldCentralH_TM( - H.xy, H.xz, H.yx, H.yz, H.zx, H.zy, surface_vertices[0][vindex], - surface_vertices[1][vindex], surface_vertices[2][vindex], &Hx, &Hy, &Hz); + CellCoordinate current_cell(surface_vertices[0][vindex], surface_vertices[1][vindex], + surface_vertices[2][vindex]); + if (params.dimension == Dimension::THREE) { + // these should adapt to use 2/1D interpolation depending on whether the y-direction is available (hence no J_tot check) + Hx = H.interpolate_to_centre_of(AxialDirection::X, current_cell); + Hy = H.interpolate_to_centre_of(AxialDirection::Y, current_cell); + Hz = H.interpolate_to_centre_of(AxialDirection::Z, current_cell); + if (J_tot != 0) { + Ex = E.interpolate_to_centre_of(AxialDirection::X, current_cell); + Ey = E.interpolate_to_centre_of(AxialDirection::Y, current_cell); + Ez = E.interpolate_to_centre_of(AxialDirection::Z, current_cell); + } else { + Ex = E.interpolate_to_centre_of(AxialDirection::X, current_cell); + Ey = E.yx[current_cell] + E.yz[current_cell]; + Ez = E.interpolate_to_centre_of(AxialDirection::Z, current_cell); + } + } else if (params.dimension == Dimension::TRANSVERSE_ELECTRIC) { + Ex = E.interpolate_to_centre_of(AxialDirection::X, current_cell); + Ey = E.interpolate_to_centre_of(AxialDirection::Y, current_cell); + Ez = 0.; + Hx = 0.; + Hy = 0.; + Hz = H.interpolate_to_centre_of(AxialDirection::Z, current_cell); + } else { + Ex = 0.; + Ey = 0.; + Ez = E.interpolate_to_centre_of(AxialDirection::Z, current_cell); + Hx = H.interpolate_to_centre_of(AxialDirection::X, current_cell); + Hy = H.interpolate_to_centre_of(AxialDirection::Y, current_cell); + Hz = 0.; + } // fprintf(stderr,"2nd interp donezn"); /*Ex and Hx*/ @@ -4960,7 +4947,7 @@ void extractPhasorsVertices(double **EHr, double **EHi, ElectricSplitField &E, M ComplexAmplitudeSample &campssample, int n, double omega, double dt, int Nt, int dimension, int J_tot, int intmethod) { - int vindex, i, j, k; + int vindex; double Ex, Ey, Ez, Hx, Hy, Hz; complex cphaseTermE, cphaseTermH; @@ -4972,48 +4959,43 @@ void extractPhasorsVertices(double **EHr, double **EHi, ElectricSplitField &E, M #pragma omp parallel default(none) \ shared(E, H, EHr, EHi, campssample) \ - private(Ex, Ey, Ez, Hx, Hy, Hz, vindex, i, j, k) \ + private(Ex, Ey, Ez, Hx, Hy, Hz, vindex) \ firstprivate(cphaseTermH, cphaseTermE, dimension, J_tot, intmethod) { #pragma omp for for (vindex = 0; vindex < campssample.n_vertices(); vindex++) { // loop over every vertex - - i = campssample.vertices[0][vindex]; - j = campssample.vertices[1][vindex]; - k = campssample.vertices[2][vindex]; - - if (dimension == THREE) - if (J_tot == 0) { - interpolateTimeDomainFieldCentralE_2Dy(E.xy, E.xz, E.yx, E.yz, E.zx, E.zy, i, j, k, &Ex, - &Ey, &Ez); - } else if (intmethod == 1) - interpolateTimeDomainFieldCentralE(E.xy, E.xz, E.yx, E.yz, E.zx, E.zy, i, j, k, &Ex, &Ey, - &Ez); - else - interpolateTimeDomainFieldCentralEBandLimited(E.xy, E.xz, E.yx, E.yz, E.zx, E.zy, i, j, k, - &Ex, &Ey, &Ez); - else if (dimension == TE) - interpolateTimeDomainFieldCentralE_TE(E.xy, E.xz, E.yx, E.yz, E.zx, E.zy, i, j, k, &Ex, &Ey, - &Ez); - else - interpolateTimeDomainFieldCentralE_TM(E.xy, E.xz, E.yx, E.yz, E.zx, E.zy, i, j, k, &Ex, &Ey, - &Ez); - if (dimension == THREE) - if (J_tot == 0) { - interpolateTimeDomainFieldCentralH_2Dy(H.xy, H.xz, H.yx, H.yz, H.zx, H.zy, i, j, k, &Hx, - &Hy, &Hz); - } else if (intmethod == 1) - interpolateTimeDomainFieldCentralH(H.xy, H.xz, H.yx, H.yz, H.zx, H.zy, i, j, k, &Hx, &Hy, - &Hz); - else - interpolateTimeDomainFieldCentralHBandLimited(H.xy, H.xz, H.yx, H.yz, H.zx, H.zy, i, j, k, - &Hx, &Hy, &Hz); - else if (dimension == TE) - interpolateTimeDomainFieldCentralH_TE(H.xy, H.xz, H.yx, H.yz, H.zx, H.zy, i, j, k, &Hx, &Hy, - &Hz); - else - interpolateTimeDomainFieldCentralH_TM(H.xy, H.xz, H.yx, H.yz, H.zx, H.zy, i, j, k, &Hx, &Hy, - &Hz); + CellCoordinate current_cell(campssample.vertices[0][vindex], campssample.vertices[1][vindex], + campssample.vertices[2][vindex]); + + if (dimension == Dimension::THREE) { + // these should adapt to use 2/1D interpolation depending on whether the y-direction is available (hence no J_tot check) + Hx = H.interpolate_to_centre_of(AxialDirection::X, current_cell); + Hy = H.interpolate_to_centre_of(AxialDirection::Y, current_cell); + Hz = H.interpolate_to_centre_of(AxialDirection::Z, current_cell); + if (J_tot != 0) { + Ex = E.interpolate_to_centre_of(AxialDirection::X, current_cell); + Ey = E.interpolate_to_centre_of(AxialDirection::Y, current_cell); + Ez = E.interpolate_to_centre_of(AxialDirection::Z, current_cell); + } else { + Ex = E.interpolate_to_centre_of(AxialDirection::X, current_cell); + Ey = E.yx[current_cell] + E.yz[current_cell]; + Ez = E.interpolate_to_centre_of(AxialDirection::Z, current_cell); + } + } else if (dimension == Dimension::TRANSVERSE_ELECTRIC) { + Ex = E.interpolate_to_centre_of(AxialDirection::X, current_cell); + Ey = E.interpolate_to_centre_of(AxialDirection::Y, current_cell); + Ez = 0.; + Hx = 0.; + Hy = 0.; + Hz = H.interpolate_to_centre_of(AxialDirection::Z, current_cell); + } else { + Ex = 0.; + Ey = 0.; + Ez = E.interpolate_to_centre_of(AxialDirection::Z, current_cell); + Hx = H.interpolate_to_centre_of(AxialDirection::X, current_cell); + Hy = H.interpolate_to_centre_of(AxialDirection::Y, current_cell); + Hz = 0.; + } update_EH(EHr, EHi, vindex, campssample.components.index(FieldComponents::Ex), cphaseTermE, Ex); update_EH(EHr, EHi, vindex, campssample.components.index(FieldComponents::Hx), cphaseTermH, Hx); diff --git a/tdms/src/main.cpp b/tdms/src/main.cpp index 3ef212339..f847bee57 100644 --- a/tdms/src/main.cpp +++ b/tdms/src/main.cpp @@ -49,18 +49,27 @@ int main(int nargs, char *argv[]){ } // decide which derivative method to use (PSTD or FDTD) - SolverMethod method = PseudoSpectral; // default + SolverMethod solver_method = PseudoSpectral; // default if (args.finite_difference()) - method = SolverMethod::FiniteDifference; + solver_method = SolverMethod::FiniteDifference; + + // decide whether to toggle off the band-limited interpolation methods + PreferredInterpolationMethods preferred_interpolation_methods = + PreferredInterpolationMethods::BandLimited;// default + if (args.cubic_interpolation()) { + preferred_interpolation_methods = PreferredInterpolationMethods::Cubic; + } //now run the time propagation code - execute_simulation(NOUTMATRICES_PASSED, (mxArray **)plhs, NMATRICES, (const mxArray **)matrixptrs, method); + execute_simulation(NOUTMATRICES_PASSED, (mxArray **) plhs, NMATRICES, + (const mxArray **) matrixptrs, solver_method, preferred_interpolation_methods); - if( !args.have_flag("-m") ){ //prints vertices and facets - saveoutput(plhs, matricestosave_all, (char **)outputmatrices_all, NOUTMATRICES_WRITE_ALL, args.output_filename()); - } - else{ // minimise the file size by not printing vertices and facets - saveoutput(plhs, matricestosave, (char **)outputmatrices, NOUTMATRICES_WRITE, args.output_filename()); + if (!args.have_flag("-m")) {//prints vertices and facets + saveoutput(plhs, matricestosave_all, (char **) outputmatrices_all, NOUTMATRICES_WRITE_ALL, + args.output_filename()); + } else {// minimise the file size by not printing vertices and facets + saveoutput(plhs, matricestosave, (char **) outputmatrices, NOUTMATRICES_WRITE, + args.output_filename()); } return 0; diff --git a/tdms/src/simulation_parameters.cpp b/tdms/src/simulation_parameters.cpp index 381af78b6..46a11a5f1 100644 --- a/tdms/src/simulation_parameters.cpp +++ b/tdms/src/simulation_parameters.cpp @@ -37,9 +37,9 @@ void SimulationParameters::set_dimension(string mode_string) { if (s == "3") { dimension = Dimension::THREE; } else if (s == "TE") { - dimension = Dimension::TE; + dimension = Dimension::TRANSVERSE_ELECTRIC; } else { - dimension = Dimension::TM; + dimension = Dimension::TRANSVERSE_MAGNETIC; } } diff --git a/tdms/tests/requirements.txt b/tdms/tests/requirements.txt index 9eec052dd..e81194ca5 100644 --- a/tdms/tests/requirements.txt +++ b/tdms/tests/requirements.txt @@ -1,3 +1,5 @@ numpy pytest +pytest-check h5py +pyyaml diff --git a/tdms/tests/system/read_config.py b/tdms/tests/system/read_config.py new file mode 100644 index 000000000..33b29376c --- /dev/null +++ b/tdms/tests/system/read_config.py @@ -0,0 +1,224 @@ +######### +# +# See YAMLTestConfig.help() for details on writing config files for tests. +# +######### +from pathlib import Path + +import yaml +from utils import run_tdms + +# config files must provide this information +required_fields = ["test_id", "tests"] +# any runs must provide this information +run_required_info = ["input_file", "reference"] +# runs may also specify this information, but if they don't we use these defaults +run_optional_info = {"fdtd_solver": False, "cubic_interpolation": False} + + +def open_and_validate_config(config_file_path) -> dict: + """Open the config file at config_file_path and validate the specifications it provides. + + Config files are required to have all required_fields present. + """ + # open the config file + config_file = open(config_file_path, "r") + # read data into dictionary + config = yaml.safe_load(config_file) + # can safely close config file after reading data + config_file.close() + + for required_field in required_fields: + if required_field not in config.keys(): + raise RuntimeError(f"{required_field} not provided in {config_file_path}") + + # pass the data from the file back + return config + + +def validate_run_information(runs: dict[str, dict[str, str]]) -> None: + """Validates the run (call to tdms) specifies all required fields. + + Each entry of "runs" contains the specifications for one run of the tdms executable. + Each such "run" is a dict that is required to have a key corresponding to the entries of run_required_info. + run_optional_info contains other properties that might be specified by a key in "run", and the defaults to use if these properties are not provided in the config file. + """ + for test_id, test_info in runs.items(): + for required_info in run_required_info: + if required_info not in test_info.keys(): + raise RuntimeError(f"{required_info} not provided in {test_id}") + return + + +class TDMSSystemTest: + """Objects of this class handle a single "run" specified in a config file, producing the command-line call to tdms and then running said command.""" + + def __init__(self, name: str, information: dict) -> None: + """Given the name of a run (EG: fs_bli) and the corresponding information about it from the config file, setup the instance.""" + # associate the command to the run name + self.run_id = name + + # check that the user has not provided any unexpected information + for info in information.keys(): + if (info not in run_required_info) and ( + info not in run_optional_info.keys() + ): + # we don't expect this information! Throw an error + raise RuntimeError(f"{self.run_id} got unexpected information: {info}") + + # extract the information about this run + # all required fields have been validated above, so they exist + for required_info in run_required_info: + setattr(self, required_info, information[required_info]) + # attempt to extract optional fields, otherwise fill with default values + for optional_info, default_option in run_optional_info.items(): + if optional_info in information.keys(): + # optional information provided + setattr(self, optional_info, information[optional_info]) + else: + # use default + setattr(self, optional_info, default_option) + return + + def get_reference_file_name(self) -> str: + """Return the name of the reference file to compare the output of this call to tdms to.""" + # reference is a required field, and an error is thrown if it is not specified. So this attribute always exists. + return self.reference + + def get_output_file_name(self) -> str: + """Returns the path that the output of the tdms executable will be saved to - this will be compared to self.reference.""" + # output files generated from running tdms are named as below + # these files are to be compared to self.reference + return f"./{self.run_id}_output.mat" + + def create_tdms_call_options(self) -> list[str]: + """Returns a list of command-line options (in appropriate order) that forms the command-line call to tdms that this instance encodes. + + This output can be passed to utils.run_tdms to run the tdms executable with the appropriate options. + """ + tdms_call_options = [] + # add optional flags + if self.fdtd_solver: + tdms_call_options.append("--finite-difference") + if self.cubic_interpolation: + tdms_call_options.append("-c") + + # add input/output file call TODO: what if we want to specify a gridfile (no tests currently do this) + tdms_call_options.append(self.input_file) + # add output file to tdms call + tdms_call_options.append(self.get_output_file_name()) + + # return a list that can be passed to run_tdms with * expansion + return tdms_call_options + + def run(self) -> None: + """Run the tdms command specified by this instance.""" + # create the options for the tdms call + tdms_command_call = self.create_tdms_call_options() + # display tdms call to the user for debugging/logging purposes + print(f"Calling tdms with arguments {tdms_command_call}") + # call run_tdms with these options + run_tdms(*tdms_command_call) + return + + +class YAMLTestConfig: + """Parser class for config.yaml files contained in .zip folders containing system test data. See the help() method for more information.""" + + def __init__(self, config_file_path=Path("./config.yaml")) -> None: + """Parse the configuration file at config_file_path.""" + # get config information + config = open_and_validate_config(config_file_path) + + # read test id + self.test_id = config["test_id"] + # read test information and create test runs + self.run_list = self.generate_test_list(config["tests"]) + + return + + def generate_test_list(self, runs: dict) -> list[TDMSSystemTest]: + """Generate a list of system tests requested in the configuration file. + + A single .zip folder, and thus configuration file, specifies multiple system tests, or "runs" of the tdms executable. Having parsed the configuration file, we can now generate a list of TDMSSystemTest objects that corresponds to these "runs", with one TDMSSystemTest object per run specified. + """ + # check that all runs have the required information + validate_run_information(runs) + # unpack the test information + run_list = [] + for run_name, run_info in runs.items(): + run_list.append(TDMSSystemTest(run_name, run_info)) + # return the list of runs (system tests) + return run_list + + def help(self) -> None: + """Prints details about input file specification to the screen.""" + help_message = """ + A single .zip file of system test data may contain multiple system tests, referred to as "runs" in this codebase so that pytest doesn't try to go crazy. Each "run" consists of a single call to the tdms executable, using a particular input file and command-line options, and then a comparison to some reference data. + + So that pytest knows which input files and options match up to particular outputs, and the tdms options that are necessary to generate them, each .zip folder contains a config.yaml with this information. These config files are stored at the top-level of the .zip folder, and use the name "config.yaml" (although there is functionality to handle using a different name if someone so wishes). + + Config file syntax is: + test_id: XX + tests: + test_name_1: + input_file: input_path_1 + reference: output_path_1 + fdtd_solver: True/False + cubic_interpolation: True/False + test_name_2: + input_file: input_path_2 + reference: output_path_1 + fdtd_solver: True/False + cubic_interpolation: True/False + Indentation is used to define blocks associated to the attribute one-indentation level higher. + + test_id is a string (although integers can also be supplied and will be converted to strings) that provides a unqiue identifier for the .zip folder. The zip folders typically follow the naming convention arc_{test_id}.zip, EG arc_01.zip and arc_example_fdtd.zip. + + Members under the "tests" identifier correspond to individual runs. + + The "test_id" and "tests" fields are required, failing to provide these will result in an error being thrown. + + Each run (member of "tests") is identified by its test_name, which must not contain spaces. These typically follow a naming convention of {solver_method}_{spatial_object} (or the reverse); for example pstd_fs or cyl_fdtd. The name can be anything you like (it is not actually used by the interpreter other than for debugging readability purposes) but for this reason we advise against naming all the runs "foobar". + + Within each run are the fields input_file, reference, and (optionally) cubic_interpolation and fdtd_solver. + The latter two are bools indicating whether the corresponding tdms option (use cubic interpolation or the fdtd solver respectively) should be toggled on. If one of these fields is missing, the default value is False (use band-limited interpolation or pstd solver methods respectively). + The input_file and reference fields are compulsory, and must specifiy the path (RELATIVE to the top level of the zipped directory) to the input file to be passed to the tdms executable, and the reference output to compare the generated output to. The generated output is named automatically. + + An example config file is below: + + '''yaml + test_id: config_example + tests: + pstd_and_bli: + input_file: input_a + reference: reference_a + pstd_and_cubic: + input_file: input_a + reference: reference_b + cubic_interpolation: True + fdtd_and_bli: + input_file: input_b + reference: reference_b + fdtd_solver: True + fdtd_and_cubic: + input_file: input_c + reference: reference_c + cubic_interpolation: True + fdtd_solver: True + ''' + + This config file corresponds to 4 runs, the tdms commands being: + - tdms input_a automatically_generated_output_name_1 + - tdms -c input_a automatically_generated_output_name_2 + - tdms --finite-difference input_b automatically_generated_output_name_3 + - tdms -c --finite-difference input_c automatically_generated_output_name_4 + + And then the comparisons made being: + - automatically_generated_output_name_1 VS reference_a + - automatically_generated_output_name_2 VS reference_b + - automatically_generated_output_name_3 VS reference_b + - automatically_generated_output_name_4 VS reference_c + """ + print(help_message) + return diff --git a/tdms/tests/system/test_arc01.py b/tdms/tests/system/test_arc01.py deleted file mode 100644 index 890cf12f4..000000000 --- a/tdms/tests/system/test_arc01.py +++ /dev/null @@ -1,37 +0,0 @@ -import os -from pathlib import Path - -from utils import HDF5File, download_data, run_tdms, work_in_zipped_dir - -ZIP_PATH = Path(os.path.dirname(os.path.abspath(__file__)), "data", "arc_01.zip") - -if not ZIP_PATH.exists(): - download_data("https://zenodo.org/record/6838866/files/arc_01.zip", to=ZIP_PATH) - - -@work_in_zipped_dir(ZIP_PATH) -def test_fs(): - """ - Ensure that the free space PSTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - """ - - run_tdms("arc_01/pstd_fs_input.mat", "pstd_fs_output.mat") - - reference = HDF5File("arc_01/pstd_fs_reference_output.mat") - assert HDF5File("pstd_fs_output.mat").matches(reference) - - -@work_in_zipped_dir(ZIP_PATH) -def test_cyl(): - """ - Ensure that the cylinder PSTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - """ - - run_tdms("arc_01/pstd_cyl_input.mat", "pstd_cyl_output.mat") - - reference = HDF5File("arc_01/pstd_cyl_reference_output.mat") - assert HDF5File("pstd_cyl_output.mat").matches(reference) diff --git a/tdms/tests/system/test_arc02.py b/tdms/tests/system/test_arc02.py deleted file mode 100644 index dc713883c..000000000 --- a/tdms/tests/system/test_arc02.py +++ /dev/null @@ -1,29 +0,0 @@ -import os -from pathlib import Path - -from utils import HDF5File, download_data, run_tdms, work_in_zipped_dir - -ZIP_PATH = Path(os.path.dirname(os.path.abspath(__file__)), "data", "arc_02.zip") - -if not ZIP_PATH.exists(): - download_data("https://zenodo.org/record/6838977/files/arc_02.zip", to=ZIP_PATH) - - -@work_in_zipped_dir(ZIP_PATH) -def test_fs(): - """Ensure that the free space PSTD output matches the reference""" - - run_tdms("arc_02/pstd_fs_input.mat", "pstd_fs_output.mat") - - reference = HDF5File("arc_02/pstd_fs_reference_output.mat") - assert HDF5File("pstd_fs_output.mat").matches(reference) - - -@work_in_zipped_dir(ZIP_PATH) -def test_cyl(): - """Ensure that the cylinder PSTD output matches the reference.""" - - run_tdms("arc_02/pstd_cyl_input.mat", "pstd_cyl_output.mat") - - reference = HDF5File("arc_02/pstd_cyl_reference_output.mat") - assert HDF5File("pstd_cyl_output.mat").matches(reference) diff --git a/tdms/tests/system/test_arc03.py b/tdms/tests/system/test_arc03.py deleted file mode 100644 index be4c16317..000000000 --- a/tdms/tests/system/test_arc03.py +++ /dev/null @@ -1,29 +0,0 @@ -import os -from pathlib import Path - -from utils import HDF5File, download_data, run_tdms, work_in_zipped_dir - -ZIP_PATH = Path(os.path.dirname(os.path.abspath(__file__)), "data", "arc_03.zip") - -if not ZIP_PATH.exists(): - download_data("https://zenodo.org/record/6839280/files/arc_03.zip", to=ZIP_PATH) - - -@work_in_zipped_dir(ZIP_PATH) -def test_fs(): - """Ensure that the free space PSTD output matches the reference""" - - run_tdms("arc_03/pstd_fs_input.mat", "pstd_fs_output.mat") - - reference = HDF5File("arc_03/pstd_fs_reference_output.mat") - assert HDF5File("pstd_fs_output.mat").matches(reference) - - -@work_in_zipped_dir(ZIP_PATH) -def test_sph(): - """Ensure that the spherical PSTD output matches the reference.""" - - run_tdms("arc_03/pstd_sph_input.mat", "pstd_sph_output.mat") - - reference = HDF5File("arc_03/pstd_sph_reference_output.mat") - assert HDF5File("pstd_sph_output.mat").matches(reference) diff --git a/tdms/tests/system/test_arc08.py b/tdms/tests/system/test_arc08.py deleted file mode 100644 index 391c86599..000000000 --- a/tdms/tests/system/test_arc08.py +++ /dev/null @@ -1,33 +0,0 @@ -import os -from pathlib import Path - -from utils import HDF5File, download_data, run_tdms, work_in_zipped_dir - -ZIP_PATH = Path(os.path.dirname(os.path.abspath(__file__)), "data", "arc_08.zip") - -if not ZIP_PATH.exists(): - download_data("https://zenodo.org/record/7147020/files/arc_08.zip", to=ZIP_PATH) - - -@work_in_zipped_dir(ZIP_PATH) -def test_fs(): - """Ensure that the free space stady-state PSTD output matches the reference""" - - run_tdms("arc_08/in/pstd_fs_steady_fast.mat", "pstd_fs_steady_output_fast.mat") - - reference = HDF5File("arc_08/out/pstd_fs_steady_fast.mat") - assert HDF5File("pstd_fs_steady_output_fast.mat").matches( - reference - ), "The free space, steady-state PSTD output doesn't match the reference." - - -@work_in_zipped_dir(ZIP_PATH) -def test_sph(): - """Ensure that the spherical steady-state PSTD output matches the reference.""" - - run_tdms("arc_08/in/pstd_sph_steady_fast.mat", "pstd_sph_steady_output_fast.mat") - - reference = HDF5File("arc_08/out/pstd_sph_steady_fast.mat") - assert HDF5File("pstd_sph_steady_output_fast.mat").matches( - reference - ), "The spherical steady-state PSTD output doesn't match the reference." diff --git a/tdms/tests/system/test_arc09.py b/tdms/tests/system/test_arc09.py deleted file mode 100644 index e172af941..000000000 --- a/tdms/tests/system/test_arc09.py +++ /dev/null @@ -1,43 +0,0 @@ -import os -from pathlib import Path - -import pytest -from utils import HDF5File, download_data, run_tdms, work_in_zipped_dir - -ZIP_PATH = Path(os.path.dirname(os.path.abspath(__file__)), "data", "arc_09.zip") - -if not ZIP_PATH.exists(): - download_data("https://zenodo.org/record/7097063/files/arc_09.zip", to=ZIP_PATH) - - -@pytest.mark.skip(reason="We know about this: likely a data problem.") -@work_in_zipped_dir(ZIP_PATH) -def test_fs(): - """ - Ensure that the free space PSTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - """ - - run_tdms("arc_09/in/pstd_fs.mat", "pstd_fs_output.mat") - - reference = HDF5File("arc_09/out/pstd_fs.mat") - assert HDF5File("pstd_fs_output.mat").matches( - reference - ), "The free space PSTD output with exdetintegral doesn't match the reference." - - -@work_in_zipped_dir(ZIP_PATH) -def test_sc(): - """ - Ensure that the (FIXME ask Peter - SC) PSTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - """ - - run_tdms("arc_09/in/pstd_sc.mat", "pstd_sc_output.mat") - - reference = HDF5File("arc_09/out/pstd_sc.mat") - assert HDF5File("pstd_sc_output.mat").matches( - reference - ), "The sc PSTD output with exdetintegral doesn't match the reference." diff --git a/tdms/tests/system/test_arc10.py b/tdms/tests/system/test_arc10.py deleted file mode 100644 index 7a4a9768e..000000000 --- a/tdms/tests/system/test_arc10.py +++ /dev/null @@ -1,41 +0,0 @@ -import os -from pathlib import Path - -from utils import HDF5File, download_data, run_tdms, work_in_zipped_dir - -ZIP_PATH = Path(os.path.dirname(os.path.abspath(__file__)), "data", "arc_10.zip") - -if not ZIP_PATH.exists(): - download_data("https://zenodo.org/record/7096040/files/arc_10.zip", to=ZIP_PATH) - - -@work_in_zipped_dir(ZIP_PATH) -def test_fs(): - """ - Ensure that the free space PSTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - """ - - run_tdms("arc_10/in/pstd_fs.mat", "pstd_fs_output.mat") - - reference = HDF5File("arc_10/out/pstd_fs.mat") - assert HDF5File("pstd_fs_output.mat").matches( - reference - ), "The free space PSTD output with exi present doesn't match the reference." - - -@work_in_zipped_dir(ZIP_PATH) -def test_cyl(): - """ - Ensure that the cylinder PSTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - """ - - run_tdms("arc_10/in/pstd_cyl.mat", "pstd_cyl_output.mat") - - reference = HDF5File("arc_10/out/pstd_cyl.mat") - assert HDF5File("pstd_cyl_output.mat").matches( - reference - ), "The spherical PSTD output with exi present doesn't match the reference." diff --git a/tdms/tests/system/test_arc12.py b/tdms/tests/system/test_arc12.py deleted file mode 100644 index b761bb2b4..000000000 --- a/tdms/tests/system/test_arc12.py +++ /dev/null @@ -1,41 +0,0 @@ -import os -from pathlib import Path - -from utils import HDF5File, download_data, run_tdms, work_in_zipped_dir - -ZIP_PATH = Path(os.path.dirname(os.path.abspath(__file__)), "data", "arc_12.zip") - -if not ZIP_PATH.exists(): - download_data("https://zenodo.org/record/7233444/files/arc_12.zip", to=ZIP_PATH) - - -@work_in_zipped_dir(ZIP_PATH) -def test_fs(): - """ - Ensure that the free space FDTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - """ - - run_tdms("--finite-difference", "arc_12/in/fdtd_fs.mat", "fdtd_fs_output.mat") - - reference = HDF5File("arc_12/out/fdtd_fs.mat") - assert HDF5File("fdtd_fs_output.mat").matches( - reference - ), "The free space FDTD output doesn't match the reference." - - -@work_in_zipped_dir(ZIP_PATH) -def test_cyl(): - """ - Ensure that the spherical FDTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - """ - - run_tdms("--finite-difference", "arc_12/in/fdtd_sph.mat", "fdtd_sph_output.mat") - - reference = HDF5File("arc_12/out/fdtd_sph.mat") - assert HDF5File("fdtd_sph_output.mat").matches( - reference - ), "The spherical FDTD output doesn't match the reference." diff --git a/tdms/tests/system/test_arc13.py b/tdms/tests/system/test_arc13.py deleted file mode 100644 index 25e2c8acc..000000000 --- a/tdms/tests/system/test_arc13.py +++ /dev/null @@ -1,41 +0,0 @@ -import os -from pathlib import Path - -from utils import HDF5File, download_data, run_tdms, work_in_zipped_dir - -ZIP_PATH = Path(os.path.dirname(os.path.abspath(__file__)), "data", "arc_13.zip") - -if not ZIP_PATH.exists(): - download_data("https://zenodo.org/record/7233118/files/arc_13.zip", to=ZIP_PATH) - - -@work_in_zipped_dir(ZIP_PATH) -def test_fs(): - """ - Ensure that the free space FDTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - """ - - run_tdms("--finite-difference", "arc_13/in/fdtd_fs.mat", "fdtd_fs_output.mat") - - reference = HDF5File("arc_13/out/fdtd_fs.mat") - assert HDF5File("fdtd_fs_output.mat").matches( - reference - ), "The free space FDTD output doesn't match the reference." - - -@work_in_zipped_dir(ZIP_PATH) -def test_cyl(): - """ - Ensure that the cylinder FDTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - """ - - run_tdms("--finite-difference", "arc_13/in/fdtd_cyl.mat", "fdtd_cyl_output.mat") - - reference = HDF5File("arc_13/out/fdtd_cyl.mat") - assert HDF5File("fdtd_cyl_output.mat").matches( - reference - ), "The spherical FDTD output doesn't match the reference." diff --git a/tdms/tests/system/test_fdtd.py b/tdms/tests/system/test_fdtd.py deleted file mode 100644 index 995038174..000000000 --- a/tdms/tests/system/test_fdtd.py +++ /dev/null @@ -1,49 +0,0 @@ -import os -from pathlib import Path - -from utils import HDF5File, download_data, run_tdms, work_in_zipped_dir - -ZIP_PATH = Path( - os.path.dirname(os.path.abspath(__file__)), "data", "arc_example_fdtd.zip" -) - -if not ZIP_PATH.exists(): - download_data( - "https://zenodo.org/record/7148198/files/arc_example_fdtd.zip", to=ZIP_PATH - ) - - -@work_in_zipped_dir(ZIP_PATH) -def test_fs(): - """ - Ensure that the free space FDTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - - Note: At the moment (due to MATLAB) this test only works if tdms is compiled with - """ - - # misleading filename here - run_tdms( - "--finite-difference", "arc_example_fdtd/in_pstd_fs.mat", "fdtd_fs_output.mat" - ) - - reference = HDF5File("arc_example_fdtd/out_fdtd_fs.mat") - assert HDF5File("fdtd_fs_output.mat").matches(reference) - - -@work_in_zipped_dir(ZIP_PATH) -def test_cyl(): - """ - Ensure that the cylinder FDTD output matches the reference. Checks - that the output .mat file (with a HDF5 format) contains tensors with - relative mean square values within numerical precision of the reference. - """ - - # misleading filename here - run_tdms( - "--finite-difference", "arc_example_fdtd/in_pstd_cyl.mat", "fdtd_cyl_output.mat" - ) - - reference = HDF5File("arc_example_fdtd/out_fdtd_cyl.mat") - assert HDF5File("fdtd_cyl_output.mat").matches(reference) diff --git a/tdms/tests/system/test_system.py b/tdms/tests/system/test_system.py new file mode 100644 index 000000000..c9a2fcfe9 --- /dev/null +++ b/tdms/tests/system/test_system.py @@ -0,0 +1,78 @@ +import os +from pathlib import Path + +import pytest +from pytest_check import check +from read_config import YAMLTestConfig +from utils import HDF5File, download_data, work_in_zipped_dir + +# Dataset is stored at https://zenodo.org/record/7440616/ +# ccaegra@ucl.ac.uk (William Graham, @willGraham01) has access. +ZENODO_URL = "https://zenodo.org/record/7440616/files" +# directory in which to store the downloaded zip files +ZIP_DESTINATION = Path(os.path.dirname(os.path.abspath(__file__)), "data") + +# all test cases and where to aquire their data +TEST_URLS = { + "01": ZENODO_URL + "/arc_01.zip", + "02": ZENODO_URL + "/arc_02.zip", + "03": ZENODO_URL + "/arc_03.zip", + "08": ZENODO_URL + "/arc_08.zip", + "09": ZENODO_URL + "/arc_09.zip", + "10": ZENODO_URL + "/arc_10.zip", + "12": ZENODO_URL + "/arc_12.zip", + "13": ZENODO_URL + "/arc_13.zip", + "example_fdtd": ZENODO_URL + "/arc_example_fdtd.zip", +} + + +def workflow() -> None: + """This is the workflow for a single zip file containing system tests: + - Extract all tests by reading the config file + - Run each system test specified in the config file + - Report failures / successes accordingly. + pytest-check is used to prevent automatic failure at the first test if the .zip file contains multiple tests. + """ + config_for_tests = YAMLTestConfig() + system_tests_in_this_folder = config_for_tests.run_list + + for system_test in system_tests_in_this_folder: + # run test + system_test.run() + # load reference file + reference_file = HDF5File(system_test.get_reference_file_name()) + # load output file + output_file = HDF5File(system_test.get_output_file_name()) + # assert file contents match + # but continue to run the other tests in this zip file, even if we fail + with check: + test_passed, comparison_information = output_file.matches( + reference_file, return_message=True + ) + assert ( + test_passed + ), f"In {config_for_tests.test_id} -> {system_test.run_id}: {comparison_information}" + return + + +# parameterise with each key, value in TEST_URLS (IE each test_id and zenodo url) +@pytest.mark.parametrize("test_id, url", list(TEST_URLS.items())) +def test_system(test_id, url) -> None: + """Fetches the system test data (from url) for test test_id and then calls workflow() on the downloaded data. + + Download is skipped if test data is already present in a cached folder - manually delete this cache to force a redownload of the data. + """ + print(f"\nRunning {test_id}", end=" | ") + + # the data should be at this location + ZIP_PATH = ZIP_DESTINATION / f"arc_{test_id}.zip" + # download data if not present + if not ZIP_PATH.exists(): + print(f"Downloading from {url}") + download_data(url, to=ZIP_PATH) + else: + print(f"Using cache in {ZIP_PATH}") + + # run the workflow in this zip folder + work_in_zipped_dir(ZIP_PATH)(workflow)() + return diff --git a/tdms/tests/system/utils.py b/tdms/tests/system/utils.py index c80276acc..d37e6cb21 100644 --- a/tdms/tests/system/utils.py +++ b/tdms/tests/system/utils.py @@ -75,7 +75,7 @@ def tuple_to_complex(dataset: h5py.Dataset) -> np.ndarray: return array.reshape(shape) - def matches(self, other: "HDF5File", rtol=1e-10) -> bool: + def matches(self, other: "HDF5File", rtol=1e-10, return_message=False) -> bool: """ Does this file match another? All arrays must be within a rtol to the other, where rtol is the relative difference between the two @@ -90,21 +90,26 @@ def matches(self, other: "HDF5File", rtol=1e-10) -> bool: other_value = other[key] if value.shape != other_value.shape: - print( - f"Array shape in {key} was not the same. " - f"{value.shape} ≠ {other_value.shape}" - ) - return False # Shapes did not match + # Shapes did not match + info_message = f"Array shape in {key} was not the same. {value.shape} ≠ {other_value.shape}" + if return_message: + return False, info_message + else: + return False r_ms_diff = relative_mean_squared_difference(value, other_value) if r_ms_diff > rtol: - print( - f"{key} was not within {rtol} to the reference. " - f"relative MSD = {r_ms_diff:.8f})" - ) - return False - - return True + # rms difference was too great + info_message = f"{key} was not within {rtol} to the reference. Relative MSD = {r_ms_diff:.8f}" + if return_message: + return False, info_message + else: + print(info_message) + return False + if return_message: + return True, "No differences detected" + else: + return True def relative_mean_squared_difference(a: np.ndarray, b: np.ndarray) -> float: diff --git a/tdms/tests/unit/field_tests/test_Field_interpolation.cpp b/tdms/tests/unit/field_tests/test_Field_interpolation.cpp index d7cc7efab..4afee58fc 100644 --- a/tdms/tests/unit/field_tests/test_Field_interpolation.cpp +++ b/tdms/tests/unit/field_tests/test_Field_interpolation.cpp @@ -10,6 +10,7 @@ #include #include +#include "cell_coordinate.h" #include "globals.h" #include "unit_test_utils.h" @@ -61,38 +62,27 @@ TEST_CASE("E-field interpolation check") { Nz = round(extent_z / cellDims[2]); // setup the "split" E-field components - ElectricSplitField E_split(Nx, Ny, Nz); - E_split.allocate(); + ElectricSplitField E_split(Nx-1, Ny-1, Nz-1); + E_split.allocate(); // alocates Nx, Ny, Nz memory space here + E_split.I_tot++; E_split.J_tot++; E_split.K_tot++; // correct the "number of datapoints" variable for these fields // setup for non-split field components ElectricField E(Nx, Ny, Nz); E.allocate(); - // storage for pointwise errors - 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); - - // compute the exact field and the "split field" components - // indices range from (0,0,0) to (Nx, Ny, Nz) INCLUSIVE! - for (int ii = 0; ii <= Nx; ii++) { - for (int jj = 0; jj <= Ny; jj++) { - for (int kk = 0; kk <= Nz; kk++) { - /* The coordinates of the grids are different, which throws a spanner in the works. - Recall that the "x-grid" is at the position of the (i,j,k)-th Ex field sample. - x-grid coordinates: (i,j,k) = (x_lower + i*cellDims[0], y_lower + (j-0.5)*cellDims[1], z_lower + (k-0.5)*cellDims[2]) - y-grid coordinates: (i,j,k) = (x_lower + (i-0.5)*cellDims[0], y_lower + j*cellDims[1], z_lower + (k-0.5)*cellDims[2]) - z-grid coordinates: (i,j,k) = (x_lower + (i-0.5)*cellDims[0], y_lower + (j-0.5)*cellDims[1], z_lower + k*cellDims[2]) - The components only depend on one of the coordinate directions, so we don't need to compute all 9 values. - */ + + /* Compute the exact field and the "split field" components + Ex[k][j][i] is the value of the field at position (x_lower,y_lower,z_lower) + (i,j+0.5,k+0.5)*cellDims + Ey[k][j][i] is the value of the field at position (x_lower,y_lower,z_lower) + (i+0.5,j,k+0.5)*cellDims + Ez[k][j][i] is the value of the field at position (x_lower,y_lower,z_lower) + (i+0.5,j+0.5,k)*cellDims + */ + for (int ii = 0; ii < Nx; ii++) { + for (int jj = 0; jj < Ny; jj++) { + for (int kk = 0; kk < Nz; kk++) { double x_comp_value = - field_component(y_lower + (((double) jj) + 0.5) * cellDims[1]);// Ex depends on y + field_component(y_lower + ((double) jj + 0.5) * cellDims[1]);// Ex depends on y double y_comp_value = - field_component(z_lower + (((double) kk) + 0.5) * cellDims[2]);// Ey depends on z + field_component(z_lower + ((double) kk + 0.5) * cellDims[2]);// Ey depends on z double z_comp_value = - field_component(x_lower + (((double) ii) + 0.5) * cellDims[0]);// Ez depends on x + field_component(x_lower + ((double) ii + 0.5) * cellDims[0]);// Ez depends on x // assign to "freq domain" ElectricField E.real.x[kk][jj][ii] = x_comp_value; E.imag.x[kk][jj][ii] = 0.; @@ -110,52 +100,63 @@ TEST_CASE("E-field interpolation check") { } } } - // the fields are now assigned, we next need to test the interpolation itself - // we perform Nx-1, Ny-1, Nz-1 interpolations here (for Ex, Ey, Ez respectively) so with some logic checks we can do everything in one for loop - // however, cell "1" has centre 0.5*cellDims, yet we are storing the true field value and interpolated values at index "0", so there are some index-offsets here + + /* In each axis, we will now interpolate to positions 1 through N{x,y,z}-1. Recall this means to the midpoint of the datapoints indexed by 0 and 1, then 1 and 2, ..., N{x,y,z}-3 and N{x,y,z}-2, and finally N{x,y,z}-2 and N{x,y,z}-1. + Ex_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 + Ey_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 + 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); + // now interpolate + // note that we aren't interpolating to position 0 (before 1st point) or N{x,y,z} (after last point) for (int ii = 0; ii < Nx; ii++) { for (int jj = 0; jj < Ny; jj++) { for (int kk = 0; kk < Nz; kk++) { // coordinates of the cell that we are interested in - double cell_centre[3]; - cell_centre[0] = x_lower + ((double) ii + 0.5) * cellDims[0]; - cell_centre[1] = y_lower + ((double) jj + 0.5) * cellDims[1]; - cell_centre[2] = z_lower + ((double) kk + 0.5) * cellDims[2]; - // Ex interpolation only goes up to Nx-1 - if (ii != Nx - 1) { - // exact field value - double Ex_exact = field_component(cell_centre[1]); - // split field interpolation + double x_eval_position = y_lower + ((double) jj + 0.5) * cellDims[1]; + double y_eval_position = z_lower + ((double) kk + 0.5) * cellDims[2]; + double z_eval_position = x_lower + ((double) ii + 0.5) * cellDims[0]; + // current cell index + CellCoordinate current_cell(ii, jj, kk); + + // Ex interpolation + if (ii!=0) { + double Ex_exact = field_component(x_eval_position);// Ex depends on y double Ex_split_interp = - E_split.interpolate_to_centre_of(AxialDirection::X, ii, jj, kk); - // freq domain field interpolation - take real part since using entirely real field - double Ex_interp = E.interpolate_to_centre_of(AxialDirection::X, ii, jj, kk).real(); - // compute the errors - Ex_error[kk][jj][ii] = Ex_interp - Ex_exact; - Ex_split_error[kk][jj][ii] = Ex_split_interp - Ex_exact; + E_split.interpolate_to_centre_of(AxialDirection::X, current_cell); + 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; } - // Ey interpolation only goes up to Ny-1 - if (jj != Ny - 1) { - double Ey_exact = field_component(cell_centre[2]); + + // Ey interpolation + if (jj!=0) { + double Ey_exact = field_component(y_eval_position);// Ey depends on z double Ey_split_interp = - E_split.interpolate_to_centre_of(AxialDirection::Y, ii, jj, kk); - double Ey_interp = E.interpolate_to_centre_of(AxialDirection::Y, ii, jj, kk).real(); - Ey_error[kk][jj][ii] = Ey_interp - Ey_exact; - Ey_split_error[kk][jj][ii] = Ey_split_interp - Ey_exact; + E_split.interpolate_to_centre_of(AxialDirection::Y, current_cell); + 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; } - // Ez interpolation only goes up to Nz-1 - if (kk != Nz - 1) { - double Ez_exact = field_component(cell_centre[0]); + + // Ez interpolation + if (kk!=0) { + double Ez_exact = field_component(z_eval_position);// Ez depends on x double Ez_split_interp = - E_split.interpolate_to_centre_of(AxialDirection::Z, ii, jj, kk); - double Ez_interp = E.interpolate_to_centre_of(AxialDirection::Z, ii, jj, kk).real(); - Ez_error[kk][jj][ii] = Ez_interp - Ez_exact; - Ez_split_error[kk][jj][ii] = Ez_split_interp - Ez_exact; + E_split.interpolate_to_centre_of(AxialDirection::Z, current_cell); + 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; } } } } - // compute error-matrix Frobenius norms double Ex_fro_err = Ex_error.frobenius(), Ey_fro_err = Ey_error.frobenius(), Ez_fro_err = Ez_error.frobenius(), Ex_split_fro_err = Ex_split_error.frobenius(), @@ -273,38 +274,29 @@ TEST_CASE("H-field interpolation check") { Nz = round(extent_z / cellDims[2]); // setup the "split" H-field components - MagneticSplitField H_split(Nx, Ny, Nz); - H_split.allocate(); + MagneticSplitField H_split(Nx - 1, Ny - 1, Nz - 1); + H_split.allocate();// alocates Nx, Ny, Nz memory space here // setup the non-split field components + H_split.I_tot++; + H_split.J_tot++; + H_split.K_tot++;// correct the "number of datapoints" variable for these fields MagneticField H(Nx, Ny, Nz); H.allocate(); - // storage for pointwise errors (-1 since we don't interpolate to cell 0's centre) - Tensor3D Hx_error, Hx_split_error, Hy_error, Hy_split_error, Hz_error, Hz_split_error; - Hx_error.allocate(Nz - 1, Ny - 1, Nx); - Hx_split_error.allocate(Nz - 1, Ny - 1, Nx); - Hy_error.allocate(Nz - 1, Ny, Nx - 1); - Hy_split_error.allocate(Nz - 1, Ny, Nx - 1); - Hz_error.allocate(Nz, Ny - 1, Nx - 1); - Hz_split_error.allocate(Nz, Ny - 1, Nx - 1); - // compute the exact field and the "split field" components - // indices range from (0,0,0) to (Nx, Ny, Nz) INCLUSIVE! - for (int ii = 0; ii <= Nx; ii++) { - for (int jj = 0; jj <= Ny; jj++) { - for (int kk = 0; kk <= Nz; kk++) { - /* The coordinates of the grids are different, which throws a spanner in the works. - Recall that the "x-grid" is at the position of the (i,j,k)-th Hx field sample. - x-grid coordinates: (i,j,k) = (x_lower + (i-0.5)*cellDims[0], y_lower + j*cellDims[1], z_lower + k*cellDims[2]) - y-grid coordinates: (i,j,k) = (x_lower + i*cellDims[0], y_lower + (j-0.5)*cellDims[1], z_lower + k*cellDims[2]) - z-grid coordinates: (i,j,k) = (x_lower + i*cellDims[0], y_lower + j*cellDims[1], z_lower + (k-0.5)*cellDims[2]) - The components only depend on one of the coordinate directions, so we don't need to compute all 9 values. - */ + /* Compute the exact field and the "split field" components + Hx_exact[k-1][j-1][i] is the value of the field at position (x_lower,y_lower,z_lower) + (i+0.5,j,k)*cellDims + Hy_exact[k-1][j][i-1] is the value of the field at position (x_lower,y_lower,z_lower) + (i,j+0.5,k)*cellDims + Hz_exact[k][j-1][i-1] is the value of the field at position (x_lower,y_lower,z_lower) + (i,j,k+0.5)*cellDims + */ + for (int ii = 0; ii < Nx; ii++) { + for (int jj = 0; jj < Ny; jj++) { + for (int kk = 0; kk < Nz; kk++) { double x_comp_value = - field_component(y_lower + ((double) jj) * cellDims[1]);// Hx depends on y + field_component(y_lower + ((double) jj + 1.) * cellDims[1]);// Hx depends on y double y_comp_value = - field_component(z_lower + ((double) kk) * cellDims[2]);// Hy depends on z + field_component(z_lower + ((double) kk + 1.) * cellDims[2]);// Hy depends on z double z_comp_value = - field_component(x_lower + ((double) ii) * cellDims[0]);// Hz depends on x + field_component(x_lower + ((double) ii + 1.) * cellDims[0]);// Hz depends on x // assign to "freq domain" ElectricField H.imag.x[kk][jj][ii] = x_comp_value; H.real.x[kk][jj][ii] = 0.; @@ -322,47 +314,60 @@ TEST_CASE("H-field interpolation check") { } } } - // the fields are now assigned, we next need to test the interpolation itself + + /* In each axis, we will now interpolate to positions 1 through N{x,y,z}-1. Recall this means to the midpoint of the datapoints indexed by 0 and 1, then 1 and 2, ..., N{x,y,z}-3 and N{x,y,z}-2, and finally N{x,y,z}-2 and N{x,y,z}-1. + Ex_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 + Ey_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 + 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 Hx_error, Hx_split_error, Hy_error, Hy_split_error, Hz_error, Hz_split_error; + Hx_error.allocate(Nz - 1, Ny - 1, Nx); + Hx_split_error.allocate(Nz - 1, Ny - 1, Nx); + Hy_error.allocate(Nz - 1, Ny, Nx - 1); + Hy_split_error.allocate(Nz - 1, Ny, Nx - 1); + Hz_error.allocate(Nz, Ny - 1, Nx - 1); + Hz_split_error.allocate(Nz, Ny - 1, Nx - 1); + + // now interpolate + // note that we aren't interpolating to position 0 (before 1st point) or N{x,y,z} (after last point) for (int ii = 0; ii < Nx; ii++) { for (int jj = 0; jj < Ny; jj++) { for (int kk = 0; kk < Nz; kk++) { // coordinates of the cell that we are interested in - double cell_centre[3]; - cell_centre[0] = x_lower + ((double) ii + 0.5) * cellDims[0]; - cell_centre[1] = y_lower + ((double) jj + 0.5) * cellDims[1]; - cell_centre[2] = z_lower + ((double) kk + 0.5) * cellDims[2]; - // Hx interpolation only goes up to Ny-1, Nz-1 - if ((jj != Ny - 1) && (kk != Nz - 1)) { - // exact field value - double Hx_exact = field_component(cell_centre[1]); - // split field interpolation + double x_eval_position = y_lower + ((double) jj + 0.5) * cellDims[1]; + double y_eval_position = z_lower + ((double) kk + 0.5) * cellDims[2]; + double z_eval_position = x_lower + ((double) ii + 0.5) * cellDims[0]; + // current cell index + CellCoordinate current_cell(ii, jj, kk); + + // Hx interpolation + if (jj != 0 && kk != 0) { + double Hx_exact = field_component(x_eval_position);// Hx depends on y double Hx_split_interp = - H_split.interpolate_to_centre_of(AxialDirection::X, ii, jj, kk); - // freq domain field interpolation - take real part since using entirely real field - double Hx_interp = - H.interpolate_to_centre_of(AxialDirection::X, ii, jj, kk).imag(); - // compute the errors - Hx_error[kk][jj][ii] = Hx_interp - Hx_exact; - Hx_split_error[kk][jj][ii] = Hx_split_interp - Hx_exact; + H_split.interpolate_to_centre_of(AxialDirection::X, current_cell); + 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; } - // Hy interpolation only goes up to Nx-1,Nz-1 - if ((ii != Nx - 1) && (kk != Nz - 1)) { - double Hy_exact = field_component(cell_centre[2]); + + // Hy interpolation + if (ii != 0 && kk != 0) { + double Hy_exact = field_component(y_eval_position);// Hy depends on z double Hy_split_interp = - H_split.interpolate_to_centre_of(AxialDirection::Y, ii, jj, kk); - double Hy_interp = - H.interpolate_to_centre_of(AxialDirection::Y, ii, jj, kk).imag(); - Hy_error[kk][jj][ii] = Hy_interp - Hy_exact; - Hy_split_error[kk][jj][ii] = Hy_split_interp - Hy_exact; + H_split.interpolate_to_centre_of(AxialDirection::Y, current_cell); + 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; } - // Hz interpolation only goes up to Nx-1,Ny-1 - if ((ii != Nx - 1) && (jj != Ny - 1)) { - double Hz_exact = field_component(cell_centre[0]); + + // Hz interpolation + if (ii != 0 && jj != 0) { + double Hz_exact = field_component(z_eval_position);// Hz depends on x double Hz_split_interp = - H_split.interpolate_to_centre_of(AxialDirection::Z, ii, jj, kk); - double Hz_interp = H.interpolate_to_centre_of(AxialDirection::Z, ii, jj, kk).imag(); - Hz_error[kk][jj][ii] = Hz_interp - Hz_exact; - Hz_split_error[kk][jj][ii] = Hz_split_interp - Hz_exact; + H_split.interpolate_to_centre_of(AxialDirection::Z, current_cell); + 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; } } } diff --git a/tdms/tests/unit/test_interpolation_determination.cpp b/tdms/tests/unit/test_interpolation_determination.cpp index 3e711f8a4..45dc05567 100644 --- a/tdms/tests/unit/test_interpolation_determination.cpp +++ b/tdms/tests/unit/test_interpolation_determination.cpp @@ -10,106 +10,99 @@ using namespace std; -/** - * @brief Test whether checkInterpolationPoints throws exceptions in the correct cases. - * - * This function checks whether it is possible to use cubic interpolation over a given index range. - * - * THIS FUNCTION WILL BE DEPRECIATED UPON SWITCHING TO THE BLI FRAMEWORK - */ -TEST_CASE("checkInterpolationPoints: exceptions thrown") { - // setup some fake field dimensions - int I = 6, J = 7, K = 8; - - // 2D simulations have J = 0 - no width in this dimension - // checkInterpolationPoints should not throw an exception if we provide j_l = 2 > j_u = -2, as is done for a 2D simulation (J=0), and request interpolation over the maximum i,k range - CHECK_NOTHROW(checkInterpolationPoints(2, I - 2, 2, -2, 2, K - 2, I, 0, K)); - - // conversely, it should throw an error if we try to provide j_l = 0 = j_u, which one may perceive are sensible choices, in this case - CHECK_THROWS_AS(checkInterpolationPoints(2, I - 2, 2, 0, 0, K - 2, I, 0, K), runtime_error); - - /* 3D simulation checks - checkInterpolationPoints should throw errors in each of the following cases: - - i_l < 2 with text Interpolation error: i_l too small - - j_l < 2 with text Interpolation error: j_l too small - - k_l < 2 with text Interpolation error: k_l too small - - i_u > I-2 with text Interpolation error: i_u too large - - j_u > J-2 with text Interpolation error: j_u too large - - k_u > K-2 with text Interpolation error: k_u too large */ - CHECK_THROWS_AS(checkInterpolationPoints(1, I - 2, 2, J - 2, 2, K - 2, I, J, K), runtime_error); - CHECK_THROWS_AS(checkInterpolationPoints(2, I - 1, 2, J - 2, 2, K - 2, I, J, K), runtime_error); - CHECK_THROWS_AS(checkInterpolationPoints(2, I - 2, 0, J - 2, 2, K - 2, I, J, K), runtime_error); - CHECK_THROWS_AS(checkInterpolationPoints(2, I - 2, 2, J - 0, 2, K - 2, I, J, K), runtime_error); - CHECK_THROWS_AS(checkInterpolationPoints(2, I - 1, 2, J - 2, -1, K - 2, I, J, K), runtime_error); - CHECK_THROWS_AS(checkInterpolationPoints(2, I - 1, 2, J - 2, 2, K + 1, I, J, K), runtime_error); -} - -/** - * @brief Test whether checkInterpolationPoints confirms interpolation is possible for the give cases. - * - * THIS FUNCTION WILL BE DEPRECIATED UPON SWITCHING TO THE BLI FRAMEWORK - */ -TEST_CASE("checkInterpolationPoints: check valid inputs") { - // setup some fake field dimensions - int I = 6, J = 7, K = 8; - - /* checkInterpolationPoints should throw errors in each of the following cases: - - i_l < 2 with text Interpolation error: i_l too small - - j_l < 2 with text Interpolation error: j_l too small - - k_l < 2 with text Interpolation error: k_l too small - - i_u > I-2 with text Interpolation error: i_u too large - - j_u > J-2 with text Interpolation error: j_u too large - - k_u > K-2 with text Interpolation error: k_u too large - When provided input arguments that do not satisfy any of the above - even cases such as i_l > i_u - should pass without throwing an exception */ - for (int k_l = 2; k_l < K - 2; k_l++) - for (int k_u = 2; k_u < K - 2; k_u++) - for (int j_l = 2; j_l < J - 2; j_l++) - for (int j_u = 2; j_u < J - 2; j_u++) - for (int i_l = 2; i_l < I - 2; i_l++) - for (int i_u = 2; i_u < I - 2; i_u++) - CHECK_NOTHROW(checkInterpolationPoints(i_l, i_u, j_l, j_u, k_l, k_u, I, J, K)); -} - /** * @brief Test whether best_scheme correctly determines the appropriate interpolation scheme to use, given the number of Yee cells either side of cell (i,j,k) * */ TEST_CASE("best_interp_scheme: correct interpolation chosen") { - SPDLOG_INFO("===== Testing interpolation scheme selection logic ====="); - int N = 10; - // should throw out_of_range exception if interpolation is impossible (<3 Yee cells in direction) + int N = 10; + bool all_schemes_correct = true; + + // should throw out_of_range exception if interpolation is impossible (<3 Yee cells in direction) + SECTION("Too few cells to interpolate") { REQUIRE_THROWS_AS(best_scheme(1, 0), out_of_range); REQUIRE_THROWS_AS(best_scheme(2, 0), out_of_range); REQUIRE_THROWS_AS(best_scheme(2, 1), out_of_range); - // should throw out_of_range exception if Yee cell if invalid index is requested + } + + /* Suppose we have N >= 8 Yee cells in a dimension. The program should determine: + - cell_id < 0 : Interpolation impossible + - cell_id == 0 : BAND_LIMITED_CELL_ZERO + - cell_id == 1,2,3 : Use BAND_LIMITED_(0,1,2) scheme respectively + - cell_id == 4,...,N-4 : Use BAND_LIMITED_3 scheme + - cell_id == N-3,N-2,N-1,N : Use BAND_LIMITED_(4,5,6,7) scheme respectively + - cell_if == N+1 : Interpolation impossible + */ + SECTION("Band-limited interpolation allowed") { + SPDLOG_INFO("Interpolation scheme selection: band-limited allowed"); + REQUIRE_THROWS_AS(best_scheme(N, -1), out_of_range); - REQUIRE_THROWS_AS(best_scheme(N, N), out_of_range); + all_schemes_correct = + all_schemes_correct && (best_scheme(N, 0).get_priority() == BAND_LIMITED_CELL_ZERO); + all_schemes_correct = + all_schemes_correct && (best_scheme(N, 1).get_priority() == BAND_LIMITED_0); + all_schemes_correct = + all_schemes_correct && (best_scheme(N, 2).get_priority() == BAND_LIMITED_1); + all_schemes_correct = + all_schemes_correct && (best_scheme(N, 3).get_priority() == BAND_LIMITED_2); + for (int i = 4; i <= N - 4; i++) { + all_schemes_correct = + all_schemes_correct && (best_scheme(N, i).get_priority() == BAND_LIMITED_3); + } + all_schemes_correct = + all_schemes_correct && (best_scheme(N, N - 3).get_priority() == BAND_LIMITED_4); + all_schemes_correct = + all_schemes_correct && (best_scheme(N, N - 2).get_priority() == BAND_LIMITED_5); + all_schemes_correct = + all_schemes_correct && (best_scheme(N, N - 1).get_priority() == BAND_LIMITED_6); + all_schemes_correct = + all_schemes_correct && (best_scheme(N, N).get_priority() == BAND_LIMITED_7); + REQUIRE(all_schemes_correct); + REQUIRE_THROWS_AS(best_scheme(N, N + 1), out_of_range); + } - /* Suppose we have N >= 8 Yee cells in a dimension. The program should determine: - - cell_id == 0 : Interpolation impossible (checked previously) - - cell_id == 1,2,3 : Use BAND_LIMITED_(0,1,2) scheme respectively - - cell_id == 4,...,N-4 : Use BAND_LIMITED_3 scheme - - cell_id == N-3,N-2,N-1 : Use BAND_LIMITED_(4,5,6) scheme respectively - */ - N = 10; - CHECK(best_scheme(N, 0).get_priority() == BAND_LIMITED_0); - CHECK(best_scheme(N, 1).get_priority() == BAND_LIMITED_1); - CHECK(best_scheme(N, 2).get_priority() == BAND_LIMITED_2); - for (int i = 3; i <= N - 4; i++) { - CHECK(best_scheme(N, i).get_priority() == BAND_LIMITED_3); + /* Suppose we have N >= 8 Yee cells in a dimension. Then if we are restricted to the cubic interpolation functions, the program should determine: + - cell_id <= 0 : Interpolation impossible + - cell_id == 1 : CUBIC_INTERP_FIRST + - cell_id == 2,3,...,N-2 : Use CUBIC_INTERP_MIDDLE + - cell_id == N-1 : CUBIC_INTERP_LAST + - cell_if >= N : Interpolation impossible + */ + SECTION("Only cubic interpolation") { + SPDLOG_INFO("Interpolation scheme selection: band-limited disallowed"); + PreferredInterpolationMethods pim = Cubic; + REQUIRE_THROWS_AS(best_scheme(N, 0, pim), out_of_range); + all_schemes_correct = + all_schemes_correct && (best_scheme(N, 1, pim).get_priority() == CUBIC_INTERP_FIRST); + for (int i = 2; i <= N - 2; i++) { + all_schemes_correct = + all_schemes_correct && (best_scheme(N, i, pim).get_priority() == CUBIC_INTERP_MIDDLE); } - CHECK(best_scheme(N, N - 3).get_priority() == BAND_LIMITED_4); - CHECK(best_scheme(N, N - 2).get_priority() == BAND_LIMITED_5); - CHECK(best_scheme(N, N - 1).get_priority() == BAND_LIMITED_6); + all_schemes_correct = + all_schemes_correct && (best_scheme(N, N - 1, pim).get_priority() == CUBIC_INTERP_LAST); + REQUIRE(all_schemes_correct); + REQUIRE_THROWS_AS(best_scheme(N, N, pim), out_of_range); + } - /* If 4 <= N < 8 we can still fall back on cubic interpolation - - cell_id == 0 : Interpolation impossible (checked previously) - - cell_id == 1 : Use CUBIC_FIRST - - cell_id == 2,...,N-2 : Use CUBIC_MIDDLE - - cell_id == N-1 : Use CUBIC_LAST - */ + /* If 4 <= N < 8 we fall back on cubic interpolation + - cell_id == 0 : Interpolation impossible + - cell_id == 1 : Use CUBIC_FIRST + - cell_id == 2,...,N-2 : Use CUBIC_MIDDLE + - cell_id == N-1 : Use CUBIC_LAST + */ + SECTION("Interpolation scheme selection: forced cubic interpolation") { + SPDLOG_INFO(""); N = 6; - CHECK(best_scheme(N, 0).get_priority() == CUBIC_INTERP_FIRST); - for(int i=1; i<=N-2; i++) {CHECK(best_scheme(N, i).get_priority() == CUBIC_INTERP_MIDDLE);} - CHECK(best_scheme(N, N-1).get_priority() == CUBIC_INTERP_LAST); + REQUIRE_THROWS_AS(best_scheme(N, 0), out_of_range); + all_schemes_correct = + all_schemes_correct && (best_scheme(N, 1).get_priority() == CUBIC_INTERP_FIRST); + for (int i = 2; i <= N - 2; i++) { + all_schemes_correct = + all_schemes_correct && (best_scheme(N, i).get_priority() == CUBIC_INTERP_MIDDLE); + } + all_schemes_correct = + all_schemes_correct && (best_scheme(N, N - 1).get_priority() == CUBIC_INTERP_LAST); + REQUIRE(all_schemes_correct); + REQUIRE_THROWS_AS(best_scheme(N, N), out_of_range); + } } diff --git a/tdms/tests/unit/test_interpolation_functions.cpp b/tdms/tests/unit/test_interpolation_functions.cpp index 27d4566d2..e80491a99 100644 --- a/tdms/tests/unit/test_interpolation_functions.cpp +++ b/tdms/tests/unit/test_interpolation_functions.cpp @@ -46,10 +46,6 @@ TEST_CASE("Cubic interpolation: exact for polynomials of degree <= 3") { double v1 = c0, v2 = c0, v3 = c0, v4 = c0; double v12 = c0, v23 = c0, v34 = c0; - // check old interp methods - CHECK(v12 == Approx(interp2(v1, v2, v3, v4)).epsilon(tol)); - CHECK(v23 == Approx(interp1(v1, v2, v3, v4)).epsilon(tol)); - CHECK(v34 == Approx(interp3(v1, v2, v3, v4)).epsilon(tol)); // check Interp_scheme class method CHECK(v12 == Approx(CBFst.interpolate(interp_data)).epsilon(tol)); CHECK(v23 == Approx(CBMid.interpolate(interp_data)).epsilon(tol)); @@ -69,10 +65,6 @@ TEST_CASE("Cubic interpolation: exact for polynomials of degree <= 3") { v23 += c1 * (x[2] + x[1]) / 2.; v34 += c1 * (x[3] + x[2]) / 2.; - // check old interp methods - CHECK(v12 == Approx(interp2(v1, v2, v3, v4)).epsilon(tol)); - CHECK(v23 == Approx(interp1(v1, v2, v3, v4)).epsilon(tol)); - CHECK(v34 == Approx(interp3(v1, v2, v3, v4)).epsilon(tol)); // check Interp_scheme class method CHECK(v12 == Approx(CBFst.interpolate(interp_data)).epsilon(tol)); CHECK(v23 == Approx(CBMid.interpolate(interp_data)).epsilon(tol)); @@ -92,10 +84,6 @@ TEST_CASE("Cubic interpolation: exact for polynomials of degree <= 3") { v23 += c2 * (x[2] + x[1]) * (x[2] + x[1]) / 4.; v34 += c2 * (x[3] + x[2]) * (x[3] + x[2]) / 4.; - // check old interp methods - CHECK(v12 == Approx(interp2(v1, v2, v3, v4)).epsilon(tol)); - CHECK(v23 == Approx(interp1(v1, v2, v3, v4)).epsilon(tol)); - CHECK(v34 == Approx(interp3(v1, v2, v3, v4)).epsilon(tol)); // check Interp_scheme class method CHECK(v12 == Approx(CBFst.interpolate(interp_data)).epsilon(tol)); CHECK(v23 == Approx(CBMid.interpolate(interp_data)).epsilon(tol)); @@ -115,10 +103,6 @@ TEST_CASE("Cubic interpolation: exact for polynomials of degree <= 3") { v23 += c3 * (x[2] + x[1]) * (x[2] + x[1]) * (x[2] + x[1]) / 8.; v34 += c3 * (x[3] + x[2]) * (x[3] + x[2]) * (x[3] + x[2]) / 8.; - // check old interp methods - CHECK(v12 == Approx(interp2(v1, v2, v3, v4)).epsilon(tol)); - CHECK(v23 == Approx(interp1(v1, v2, v3, v4)).epsilon(tol)); - CHECK(v34 == Approx(interp3(v1, v2, v3, v4)).epsilon(tol)); // check Interp_scheme class method CHECK(v12 == Approx(CBFst.interpolate(interp_data)).epsilon(tol)); CHECK(v23 == Approx(CBMid.interpolate(interp_data)).epsilon(tol)); @@ -280,7 +264,7 @@ TEST_CASE("BLi: MATLAB benchmarking") { f_data[nSamples - 1] = analytic_function(xi[nSamples - 1]); // perform interpolation for (int i = 0; i < nSamples - 1; i++) { - InterpolationScheme scheme = best_scheme(nSamples - 1, i); + InterpolationScheme scheme = best_scheme(nSamples, i + 1); f_interp[i] = scheme.interpolate(f_data, i + 1 - scheme.number_of_datapoints_to_left); // Compare interpolated values to the true values f_errors[i] = abs(f_exact[i] - f_interp[i]); @@ -303,7 +287,7 @@ TEST_CASE("BLi: MATLAB benchmarking") { f_data[nSamples - 1] = complex_fn(xi[nSamples - 1]); // perform interpolation for (int i = 0; i < nSamples - 1; i++) { - InterpolationScheme scheme = best_scheme(nSamples - 1, i); + InterpolationScheme scheme = best_scheme(nSamples, i + 1); f_interp[i] = scheme.interpolate(f_data, i + 1 - scheme.number_of_datapoints_to_left); // Compare interpolated values to the true values f_errors[i] = abs(f_exact[i] - f_interp[i]); @@ -315,7 +299,6 @@ TEST_CASE("BLi: MATLAB benchmarking") { REQUIRE(is_close_or_better(max_error, MATLAB_max_error)); // Compare norm-errors norm_error = euclidean(f_errors, nSamples - 2); - fprintf(stderr, "%.8e vs %.8e\n", norm_error, MATLAB_norm_error); REQUIRE(is_close_or_better(norm_error, MATLAB_norm_error)); // report test results