Merge branch 'stuart-knock-wave_to_rennie'

This pull request solves #80 and #45.

Author: pausz

src/harmonic.cpp

@@ -25,7 +25,7 @@ using std::vector;
 /**
   @brief Computes the derivatives of the Harmonic Propagator function.
 
-  The Harmonic propgator equation is given by:
+  The Harmonic propagator equation is given by:
   \f{eqnarray*}{
     \frac{d^2\phi}{dt^2}&=& \gamma^2 \left(Q - \frac{2}{\gamma} \frac{d\phi}{dt} - \phi \right)
   \f}

src/solver.cpp

@@ -1,7 +1,12 @@
 /** @file solver.cpp
-  @brief A brief, one sentence description.
+  @brief Solver member function definitions.
 
-  A more detailed multiline description...
+  The Solver class pulls together Population, Propagator, and Coupling classes.
+
+  If you define your own Propagators or Couplings, then, it is in the init
+  member-function of this class that you must create appropriate entries in
+  order to make those Propagators and Couplings specifiable via neurofield's
+  configuration (.conf) files.
 
   @author Peter Drysdale, Felix Fung,
 */
@@ -23,6 +28,7 @@
 #include "output.h"        // Output; Outlet;
 #include "population.h"    // Population;
 #include "wave.h"          // Wave;
+#include "wave_legacy.h"   // WaveLegacy;
 
 // C++ standard library headers
 #include <cmath>    // std::remainder; std::sqrt;
@@ -162,6 +168,14 @@ void Solver::init( Configf& configf ) {
     } else if(ptype=="HarmonicIntegral") {
       propagators.add( new
                    HarmonicIntegral(nodes,deltat,i, *pops[cnt.pre[i]], *pops[cnt.post[i]], longside, topology));
+    } else if(ptype=="WaveLegacy") {
+      if( nodes==1 ) {
+        propagators.add( new
+                     Harmonic(nodes,deltat,i, *pops[cnt.pre[i]], *pops[cnt.post[i]], longside, topology));
+      } else {
+        propagators.add( new
+                     WaveLegacy(nodes,deltat,i, *pops[cnt.pre[i]], *pops[cnt.post[i]], longside, topology));
+      }
     } else {
       cerr<<"Invalid propagator type '"<<ptype<<"'."<<endl;
       exit(EXIT_FAILURE);

src/solver.h

@@ -1,5 +1,5 @@
 /** @file solver.h
-  @brief Defines the Solver class, a high level container for managing a simulation.
+  @brief Declares the Solver class, a high level container for managing a simulation.
 
   The Solver class pulls together Population, Propagator, and Coupling classes.
 

src/stencil.cpp

@@ -148,7 +148,6 @@ const vector<double>& Stencil::operator= ( const vector<double>& field ) {
   return field;
 }
 
-
 void Stencil::set( int node ) const {
   if( node>=0 && node<nodes ) {
     int x = node%longside;

src/stencil.h

@@ -1,7 +1,19 @@
 /** @file stencil.h
-  @brief A brief, one sentence description.
+  @brief Stencil declaration -- spatial neighbourhood use to evaluate \f[{\Delta}x\f].
 
-  A more detailed multiline description...
+  This file contains the Stencil declaration, including a definition of Stencil's
+  grid increment operator (ie, ++), and the parentheses operator (ie, ()). The
+  default neighbourhood is the orthogonal or von Neumann neighbourhood of a
+  two-dimensional five-point stencil, that is:
+  \f{matrix}{
+      & n &   \\
+    w & c & e \\
+      & s
+  \f}
+
+  #TODO: a reference to a good review paper or textbook would be better here.
+  https://en.wikipedia.org/wiki/Five-point_stencil
+  https://en.wikipedia.org/wiki/Stencil_code
 
   @author Peter Drysdale, Felix Fung,
 */
@@ -13,12 +25,6 @@
 #include <string> // std::string;
 #include <vector> // std::vector;
 
-enum Moore {
-  nw=-4, n, ne,
-  w,     c,  e,
-  sw,    s, se
-};
-
 class Stencil {
   Stencil();
   Stencil& operator=(const Stencil&);
@@ -32,27 +38,32 @@ class Stencil {
 
   mutable int ptr;
  public:
+  enum Neighbours {
+       n=-2,
+    w, c, e,
+       s
+  };
+
   Stencil( int nodes, int longside, const std::string& boundary );
   virtual ~Stencil(void);
 
   const std::vector<double>& operator= ( const std::vector<double>& field );
-  inline void operator++ (int) const; // increment Moore grid
+  inline void operator++ (int) const; // increment Neighbours grid
   void set( int node ) const; // point to node
   int get(void) const; // get ptr
 
-  inline double operator() ( Moore moore=c ) const {
-    int arr[9] = {ptr-longside-2-1, ptr-longside-2, ptr-longside-2+1,
-                  ptr-1,            ptr,            ptr+1,
-                  ptr+longside+2-1, ptr+longside+2, ptr+longside+2+1};
-    return m[arr[moore+4]];
+  inline double operator() ( Neighbours neighbour=c ) const {
+    int arr[5] = {       ptr-longside-2,
+                  ptr-1,       ptr,      ptr+1,
+                         ptr+longside+2};
+    return m[arr[neighbour+2]];
   }
 
   inline operator double (void) const {
     return (*this)(c);
   }
 };
 
-
 void Stencil::operator++ (int) const {
   ptr++;
   if( ( ptr%(longside+2)==longside+1 ) ) {

src/stencil_legacy.cpp

@@ -0,0 +1,177 @@
+/** @file stencil_legacy.cpp
+  @brief StencilLegacy, square nine-point stencil definition.
+
+  As well as a constructor and destructor, this file defines the set() and
+  get() member functions and a custom equals (=) operator.
+
+  @author Peter Drysdale, Felix Fung.
+*/
+
+// Main module header
+#include "stencil_legacy.h"  // StencilLegacy
+
+// C++ standard library headers
+#include <iostream> // std::cerr; std::endl;
+#include <string>   // std::string;
+#include <vector>   // std::vector;
+using std::cerr;
+using std::endl;
+using std::string;
+using std::vector;
+
+StencilLegacy::StencilLegacy( int nodes, int longside, const string& boundary )
+  : nodes(nodes), longside(longside), shortside(nodes/longside),
+    boundary(boundary), ptr(0) {
+  if( boundary != "Torus" && boundary != "Sphere" ) {
+    cerr<<"StencilLegacy boundary must be Torus or Sphere."<<endl;
+    exit(EXIT_FAILURE);
+  }
+  m = new double[ (shortside+2)*(longside+2) ];
+}
+
+StencilLegacy::~StencilLegacy() {
+  delete[] m;
+}
+
+const vector<double>& StencilLegacy::operator= ( const vector<double>& field ) {
+  if( boundary == "Torus" ) {
+    // copy centre square
+    const double* __restrict__ p1 = &field[0];
+    double* __restrict__ p2 = &m[3+longside];
+    for( int i=0; i<shortside; i++ ) {
+      for( int j=0; j<longside; j++ ) {
+        *p2++ = *p1++;
+      }
+      p2 += 2;
+    }
+
+    // copy top row from bottom row
+    p1 = &field[longside*(shortside-1)];
+    p2 = &m[1];
+    for( int i=0; i<longside; i++ ) {
+      *p2++ = *p1++;
+    }
+
+    // copy bottom row from top row
+    p1 = &field[0];
+    p2 = &m[(longside+2)*(shortside+1) +1];
+    for( int i=0; i<longside; i++ ) {
+      *p2++ = *p1++;
+    }
+
+    // copy left edge from right edge
+    p1 = &field[longside -1];
+    p2 = &m[longside +2];
+    for( int i=0; i<shortside; i++, p1+=longside, p2+=longside+2 ) {
+      *p2 = *p1;
+    }
+
+    // copy right edge from left edge
+    p1 = &field[0];
+    p2 = &m[2*longside +3];
+    for( int i=0; i<shortside; i++, p1+=longside, p2+=longside+2 ) {
+      *p2 = *p1;
+    }
+
+    /*// copy centre square
+    for( int j=0; j< nodes/longside; j++ )
+      for( int i=0; i<longside; i++ )
+        m[(j+1)*(longside+2)+i+1] = field[j*longside+i];
+
+    // copy right edge into left boundary
+    for( int i=0; i<=nodes/longside; i++ )
+      m[i*(longside+2)] = field[i*longside-1];
+
+    // copy left edge into right boundary
+    for( int i=0; i<=nodes/longside; i++ )
+      m[(i+2)*(longside+2)-1] = field[i*longside];
+
+    // copy bottom edge into top boundary
+    for( int i=0; i<longside; i++ )
+      m[i+1] = field[nodes-longside+i];
+
+    // copy top edge into bottom boundary
+    for( int i=0; i<longside; i++ )
+      m[(nodes/longside+1)*(longside+2)+1+i] = field[i];*/
+
+    // copy 4 corners
+    m[0] = field[nodes-1];
+    m[longside+1] = field[nodes-longside];
+    m[(longside+2)*(nodes/longside+1)] = field[longside-1];
+    m[(longside+2)*(nodes/longside+2)-1] = field[0];
+  } else if( boundary == "Sphere" ) {
+    // copy centre square
+    const double* __restrict__ p1 = &field[0];
+    double* __restrict__ p2 = &m[3+longside];
+    for( int i=0; i<shortside; i++ ) {
+      for( int j=0; j<longside; j++ ) {
+        *p2++ = *p1++;
+      }
+      p2 += 2;
+    }
+
+    // copy top row from bottom row
+    p1 = &field[longside*(shortside-1)];
+    p2 = &m[longside+1];
+    for( int i=0; i<longside; i++ ) {
+      *p2-- = *p1++;
+    }
+
+    // copy bottom row from top row
+    p1 = &field[longside-1];
+    p2 = &m[(longside+2)*(shortside+1) +1];
+    for( int i=0; i<longside; i++ ) {
+      *p2++ = *p1--;
+    }
+
+    // copy left edge from right edge
+    p1 = &field[longside -1];
+    p2 = &m[longside +2];
+    for( int i=0; i<shortside; i++, p1+=longside, p2+=longside+2 ) {
+      *p2 = *p1;
+    }
+
+    // copy right edge from left edge
+    p1 = &field[0];
+    p2 = &m[2*longside +3];
+    for( int i=0; i<shortside; i++, p1+=longside, p2+=longside+2 ) {
+      *p2 = *p1;
+    }
+
+    // copy 4 corners NEEDS TO MODIFY SOMEWHERE!!
+    m[0] = field[nodes-1];
+    m[longside+1] = field[nodes-longside];
+    m[(longside+2)*(nodes/longside+1)] = field[longside-1];
+    m[(longside+2)*(nodes/longside+2)-1] = field[0];
+  }
+
+  set(0);
+  return field;
+}
+
+void StencilLegacy::set( int node ) const {
+  if( node>=0 && node<nodes ) {
+    int x = node%longside;
+    int y = node/longside;
+    ptr = (y+1)*(longside+2) +x+1;
+  } else {
+    cerr<<"StencilLegacy node setting out of bound: "<<node<<endl;
+    exit(EXIT_FAILURE);
+  }
+}
+
+int StencilLegacy::get() const {
+  int x = ptr%(longside+2)-1;
+  int y = ptr/(longside+2)-1;
+  return x+y*longside;
+}
+
+/*double StencilLegacy::nw(void) const { return m[ ptr-longside-2-1]; }
+double StencilLegacy::n (void) const { return m[ ptr-longside-2  ]; }
+double StencilLegacy::ne(void) const { return m[ ptr-longside-2+1]; }
+double StencilLegacy:: w(void) const { return m[ ptr           -1]; }
+double StencilLegacy:: c(void) const { return m[ ptr             ]; }
+double StencilLegacy:: e(void) const { return m[ ptr           +1]; }
+double StencilLegacy::sw(void) const { return m[ ptr+longside+2-1]; }
+double StencilLegacy::s (void) const { return m[ ptr+longside+2  ]; }
+double StencilLegacy::se(void) const { return m[ ptr+longside+2+1]; }*/

src/stencil_legacy.h

@@ -0,0 +1,74 @@
+/** @file stencil_legacy.h
+  @brief StencilLegacy square nine-point stencil declaration.
+
+  The StencilLegacy stencil declaration, including definition of the Moore
+  neighbourhood and the stencil's Moore-grid increment operator (ie, ++).
+  The specific Moore neighbourhood used here implements a "square" nine-point
+  stencil.
+
+  #TODO: a reference to a good review paper or textbook would be better here.
+  https://en.wikipedia.org/wiki/Stencil_code
+
+  @author Peter Drysdale, Felix Fung.
+*/
+
+#ifndef NEUROFIELD_SRC_STENCIL_LEGACY_H
+#define NEUROFIELD_SRC_STENCIL_LEGACY_H
+
+// C++ standard library headers
+#include <string> // std::string;
+#include <vector> // std::vector;
+
+class StencilLegacy {
+  StencilLegacy();
+  StencilLegacy& operator=(const StencilLegacy&);
+  StencilLegacy(const StencilLegacy&);
+ protected:
+  int nodes;
+  int longside;
+  int shortside;
+  std::string boundary;
+  double* m;
+
+  mutable int ptr;
+ public:
+  enum Moore {
+    nw=-4, n, ne,
+    w,     c,  e,
+    sw,    s, se
+  };
+
+  StencilLegacy( int nodes, int longside, const std::string& boundary );
+  virtual ~StencilLegacy(void);
+
+  const std::vector<double>& operator= ( const std::vector<double>& field );
+  inline void operator++ (int) const; // increment Moore grid
+  void set( int node ) const; // point to node
+  int get(void) const; // get ptr
+
+  inline double operator() ( Moore moore=c ) const {
+    int arr[9] = {ptr-longside-2-1, ptr-longside-2, ptr-longside-2+1,
+                  ptr-1,            ptr,            ptr+1,
+                  ptr+longside+2-1, ptr+longside+2, ptr+longside+2+1};
+    return m[arr[moore+4]];
+  }
+
+  inline operator double (void) const {
+    return (*this)(c);
+  }
+};
+
+
+void StencilLegacy::operator++ (int) const {
+  ptr++;
+  if( ( ptr%(longside+2)==longside+1 ) ) {
+    // Move stencil from far side to near side
+    ptr += 2;
+  }
+  if( ptr == (longside+1)*(shortside+2)+1 ) {
+    // Move stencil from very end of the array to the
+    set(0);
+  }
+}
+
+#endif //NEUROFIELD_SRC_STENCIL_LEGACY_H