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 <[email protected]> * Update tdms/include/field.h Missed one Co-authored-by: Sam Cunliffe <[email protected]> * 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 <[email protected]> * 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 <[email protected]>
UCL · Dec 16, 2022 · a60f7be · a60f7be
1 parent 063f4a1
commit a60f7be
Show file tree

Hide file tree

Showing 35 changed files with 1,207 additions and 2,389 deletions.
diff --git 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
@@ -10,6 +10,7 @@
 
 #include <fftw3.h>
 
+#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
@@ -3,6 +3,7 @@
  * @author William Graham ([email protected])
  * @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
@@ -6,6 +6,7 @@
 
 #include <complex>
 
+#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<double>{
@@ -45,6 +54,9 @@ class SplitFieldComponent: public Tensor3D<double>{
   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,15 +198,14 @@ 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.; };
 };
 
 /**
  * A complex field defined over a grid. Has real and imaginary components
  * 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<double> The interpolated field value
    */
-  virtual std::complex<double> interpolate_to_centre_of(AxialDirection d, int i, int j, int k) = 0;
+  virtual std::complex<double> 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<double> *x_at_centre,
+                                                          std::complex<double> *y_at_centre,
+                                                          std::complex<double> *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<double> *x_at_centre,
+                                                          std::complex<double> *y_at_centre,
+                                                          std::complex<double> *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<double> The interpolated component value
    */
-  std::complex<double> interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override;
+  std::complex<double> 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<double> *x_at_centre,
+                                                  std::complex<double> *y_at_centre,
+                                                  std::complex<double> *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<double> *x_at_centre,
+                                                  std::complex<double> *y_at_centre,
+                                                  std::complex<double> *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<double> The interpolated component value
    */
-  std::complex<double> interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override;
+  std::complex<double> 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<double> *x_at_centre,
+                                                  std::complex<double> *y_at_centre,
+                                                  std::complex<double> *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<double> *x_at_centre,
+                                                  std::complex<double> *y_at_centre,
+                                                  std::complex<double> *z_at_centre) override;
 };
 
 /**

diff --git 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