E- and H-field interpolation tests, frequency domain interpolation

tdms/CMakeLists.txt

@@ -154,6 +154,8 @@ if (BUILD_TESTING)
 
     add_executable(tdms_tests
             tests/unit/test_fields.cpp
+            tests/unit/test_BLi_vs_cubic_interpolation.cpp
+            tests/unit/test_field_interpolation.cpp
             tests/unit/test_interpolation_determination.cpp
             tests/unit/test_interpolation_functions.cpp
             tests/unit/test_numerical_derivative.cpp

tdms/include/arrays.h

@@ -25,6 +25,13 @@ class XYZTensor3D {
       default: throw std::runtime_error("Have no element " + to_string(c));
     }
   }
+
+  /**
+   * @brief Allocates x, y, and z as (K_total+1) * (J_total+1) * (I_total+1) arrays
+   * 
+   * @param I_total,J_total,K_total Dimensions of the tensor size to set
+   */
+  void allocate(int I_total, int J_total, int K_total);
 };
 
 class XYZVectors {
@@ -222,7 +229,7 @@ class Tensor3D{
 
   void allocate(int nK, int nJ, int nI){
     n_layers = nK, n_cols = nJ, n_rows = nI;
-    tensor = (T ***)malloc(n_layers * sizeof(T *));
+    tensor = (T ***)malloc(n_layers * sizeof(T **));
 
     for(int k=0; k < n_layers; k++){
       tensor[k] = (T **)malloc(n_cols * sizeof(T *));
@@ -235,6 +242,13 @@ class Tensor3D{
     }
   };
 
+  /**
+   * @brief Computes the Frobenius norm of the tensor
+   * 
+   * fro_norm = \f$\sqrt{ \sum_{i=0}^{I_tot}\sum_{j=0}^{J_tot}\sum_{k=0}^{K_tot} |t[k][j][i]|^2 }\f$
+   */
+  double frobenius();
+
   ~Tensor3D(){
     if (tensor == nullptr) return;
     if (is_matlab_initialised){

tdms/include/field.h

@@ -10,6 +10,7 @@
 #include "mat_io.h"
 #include "simulation_parameters.h"
 #include "utils.h"
+#include "globals.h"
 
 
 /**
@@ -22,6 +23,9 @@
  *  (0, 0, 0)   (1, 0, 0)     (2, 0, 0)
  *
  *  has I_tot = 2, J_tot = 1, K_tot = 0.
+ * 
+ * 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{
 
@@ -111,6 +115,15 @@ class SplitField : public Grid{
      * @param eh_vec // TODO
      */
     void initialise_fftw_plan(int n_threads, EHVec &eh_vec);
+
+    /**
+   * @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 
+   * @return double The interpolated field value
+   */
+    virtual double interpolate_to_centre_of(AxialDirection d, int i, int j, int k) = 0;
 };
 
 class ElectricSplitField: public SplitField{
@@ -126,6 +139,15 @@ class ElectricSplitField: public SplitField{
      */
     ElectricSplitField(int I_total, int J_total, int K_total) :
             SplitField(I_total, J_total, K_total){};
+
+    /**
+   * @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 
+   * @return double The interpolated component value
+   */
+    double interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override;
 };
 
 class MagneticSplitField: public SplitField{
@@ -141,6 +163,15 @@ class MagneticSplitField: public SplitField{
      */
     MagneticSplitField(int I_total, int J_total, int K_total) :
             SplitField(I_total, J_total, K_total){};
+
+    /**
+   * @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 
+   * @return double The interpolated component value
+   */
+    double interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override;
 };
 
 class CurrentDensitySplitField: public SplitField{
@@ -156,6 +187,8 @@ 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.; };
 };
 
 /**
@@ -169,7 +202,7 @@ class Field : public Grid{
 
   std::complex<double> angular_norm = 0.;
 
-  // TODO: this is likely better as a set of complex arrays
+  // TODO: this is likely better as a set of complex arrays - use Tensor3d<std::complex<double>>
   XYZTensor3D real;
   XYZTensor3D imag;
 
@@ -184,10 +217,35 @@ class Field : public Grid{
   void normalise_volume();
 
   /**
-     * Zero all components of the real and imaginary parts of the field
+     * Default no arguments constructor
+     */
+  Field() = default;
+
+  /**
+     * Constructor of the field with a defined size in the x, y, z Cartesian
+     * dimensions
+     */
+  Field(int I_total, int J_total, int K_total);
+
+  /**
+     * Allocate the memory appropriate for all the 3D tensors associated with
+     * this split field
+     */
+  void allocate();
+
+  /**
+     * Set all the values of all components of the field to zero
      */
   void zero();
 
+  /**
+     * Allocate and set to zero all components of the field
+     */
+  void allocate_and_zero() {
+    allocate();
+    zero();
+  }
+
   /**
    * Set the phasors for this field, given a split field. Result gives field according to the
    * exp(-iwt) convention
@@ -206,6 +264,14 @@ class Field : public Grid{
   std::complex<double> phasor_norm(double f, int n, double omega, double dt, int Nt);
 
   virtual double phase(int n, double omega, double dt) = 0;
+  /**
+   * @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 
+   * @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;
 
   /**
    * Set the values of all components in this field from another, equally sized field
@@ -219,12 +285,38 @@ class ElectricField: public Field{
 
 private:
   double phase(int n, double omega, double dt) override;
+
+public:
+  ElectricField() = default;
+  ElectricField(int I_total, int J_total, int K_total) : Field(I_total, J_total, K_total){};
+
+  /**
+   * @brief Interpolates an 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 
+   * @return std::complex<double> The interpolated component value
+   */
+  std::complex<double> interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override;
 };
 
 class MagneticField: public Field{
 
 private:
   double phase(int n, double omega, double dt) override;
+
+public:
+  MagneticField() = default;
+  MagneticField(int I_total, int J_total, int K_total) : Field(I_total, J_total, K_total){};
+
+  /**
+   * @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 
+   * @return std::complex<double> The interpolated component value
+   */
+  std::complex<double> interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override;
 };
 
 /**

tdms/include/globals.h

@@ -38,6 +38,12 @@ typedef struct complex_vector
 	std::complex<double> Z;
 } complex_vector;
 
+enum AxialDirection 
+{
+	X = 'x',
+	Y = 'y',
+	Z = 'z'
+};
 
 // **********************
 //  Enumerated constants

tdms/include/interpolation_methods.h

@@ -1,12 +1,14 @@
 /**
- * @file interpolation_methods.cpp
+ * @file interpolation_methods.h
  * @author William Graham ([email protected])
  * @brief InterpScheme class methods and supporting functions
  * 
  * Non InterpScheme methods are required to preserve functionality whilst testing new schemes
  */
 #pragma once
 
+#include <complex>
+
 /*Use cubic interpolation to interpolate between the middle 2 of 4 points
  * v0    v1    v2    v3
  * o     o  x   o    o
@@ -78,7 +80,7 @@ enum scheme_value
     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]
+    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)
@@ -129,14 +131,24 @@ class InterpolationScheme {
          * @brief Executes the interpolation scheme on the data provided
          * 
          * The interpolation schemes are all of the form
-         * interpolated_value = \sum_{i=first_nonzero_coeff}^{last_nonzero_coeff} scheme_coeffs[i] * v[i],
+         * interpolated_value = \sum_{i=0}^{7} scheme_coeffs[i] * v[i],
          * so provided that the coefficients have been set correctly in construction (and the data gathered appropriately), we can run the same for loop for each interpolation scheme.
          * 
-         * @param v Sample datapoints to use in interpolation
+         * For slight speedup, the actual sum performed loops over those i such that
+         * 0 <= first_nonzero_coeff <= i <= last_nonzero_coeff <= 7.
+         * 
+         * @param v Sample datapoints to use in interpolation; v[0] should be the first of 8 values
          * @param offset [Default 0] Read buffer from v[offset] rather than v[0]
          * @return double Interpolated value
          */
-        double interpolate(const double *v, const int offset = 0) const;
+        template<typename T>
+        T interpolate(const T *v, const int offset = 0) const {
+          T interp_value = 0.;
+          for (int ind = first_nonzero_coeff; ind <= last_nonzero_coeff; ind++) {
+            interp_value += scheme_coeffs[ind] * v[ind + offset];
+          }
+          return interp_value;
+        };
 
         /**
          * @brief Determines whether another interpScheme has greater value than this one