Regenerate code

6 files changed, 87 insertions, 162 deletions

diff --git a/ML_WaveToyFO/interface.ccl b/ML_WaveToyFO/interface.ccl
index a9f7ff5..cdb6319 100644
--- a/ML_WaveToyFO/interface.ccl
+++ b/ML_WaveToyFO/interface.ccl
@@ -8,6 +8,7 @@ inherits:   Grid GenericFD Boundary
 
 USES INCLUDE: GenericFD.h
 USES INCLUDE: Symmetry.h
+USES INCLUDE: Vectors.hh
 USES INCLUDE: sbp_calc_coeffs.h
 USES INCLUDE: Boundary.h
 USES INCLUDE: loopcontrol.h
diff --git a/ML_WaveToyFO/src/Differencing.h b/ML_WaveToyFO/src/Differencing.h
index aa21060..1c7a87f 100644
--- a/ML_WaveToyFO/src/Differencing.h
+++ b/ML_WaveToyFO/src/Differencing.h
@@ -1,12 +1,12 @@
-#define PDstandardNth1(u,i,j,k) (p1o12dx*(-8*(u)[index+di*(-1)+dj*(0)+dk*(0)] + 8*(u)[index+di*(1)+dj*(0)+dk*(0)] + (u)[index+di*(-2)+dj*(0)+dk*(0)] - (u)[index+di*(2)+dj*(0)+dk*(0)]))
-#define PDstandardNth2(u,i,j,k) (p1o12dy*(-8*(u)[index+di*(0)+dj*(-1)+dk*(0)] + 8*(u)[index+di*(0)+dj*(1)+dk*(0)] + (u)[index+di*(0)+dj*(-2)+dk*(0)] - (u)[index+di*(0)+dj*(2)+dk*(0)]))
-#define PDstandardNth3(u,i,j,k) (p1o12dz*(-8*(u)[index+di*(0)+dj*(0)+dk*(-1)] + 8*(u)[index+di*(0)+dj*(0)+dk*(1)] + (u)[index+di*(0)+dj*(0)+dk*(-2)] - (u)[index+di*(0)+dj*(0)+dk*(2)]))
-#define PDstandardNth11(u,i,j,k) (pm1o12dx2*(30*(u)[index+di*(0)+dj*(0)+dk*(0)] - 16*((u)[index+di*(-1)+dj*(0)+dk*(0)] + (u)[index+di*(1)+dj*(0)+dk*(0)]) + (u)[index+di*(-2)+dj*(0)+dk*(0)] + (u)[index+di*(2)+dj*(0)+dk*(0)]))
-#define PDstandardNth22(u,i,j,k) (pm1o12dy2*(30*(u)[index+di*(0)+dj*(0)+dk*(0)] - 16*((u)[index+di*(0)+dj*(-1)+dk*(0)] + (u)[index+di*(0)+dj*(1)+dk*(0)]) + (u)[index+di*(0)+dj*(-2)+dk*(0)] + (u)[index+di*(0)+dj*(2)+dk*(0)]))
-#define PDstandardNth33(u,i,j,k) (pm1o12dz2*(30*(u)[index+di*(0)+dj*(0)+dk*(0)] - 16*((u)[index+di*(0)+dj*(0)+dk*(-1)] + (u)[index+di*(0)+dj*(0)+dk*(1)]) + (u)[index+di*(0)+dj*(0)+dk*(-2)] + (u)[index+di*(0)+dj*(0)+dk*(2)]))
-#define PDstandardNth12(u,i,j,k) (p1o144dxdy*(-64*((u)[index+di*(-1)+dj*(1)+dk*(0)] + (u)[index+di*(1)+dj*(-1)+dk*(0)]) + 64*((u)[index+di*(-1)+dj*(-1)+dk*(0)] + (u)[index+di*(1)+dj*(1)+dk*(0)]) + 8*((u)[index+di*(-1)+dj*(2)+dk*(0)] + (u)[index+di*(1)+dj*(-2)+dk*(0)] + (u)[index+di*(-2)+dj*(1)+dk*(0)] + (u)[index+di*(2)+dj*(-1)+dk*(0)]) - 8*((u)[index+di*(-1)+dj*(-2)+dk*(0)] + (u)[index+di*(1)+dj*(2)+dk*(0)] + (u)[index+di*(-2)+dj*(-1)+dk*(0)] + (u)[index+di*(2)+dj*(1)+dk*(0)]) + (u)[index+di*(-2)+dj*(-2)+dk*(0)] - (u)[index+di*(-2)+dj*(2)+dk*(0)] - (u)[index+di*(2)+dj*(-2)+dk*(0)] + (u)[index+di*(2)+dj*(2)+dk*(0)]))
-#define PDstandardNth13(u,i,j,k) (p1o144dxdz*(-64*((u)[index+di*(-1)+dj*(0)+dk*(1)] + (u)[index+di*(1)+dj*(0)+dk*(-1)]) + 64*((u)[index+di*(-1)+dj*(0)+dk*(-1)] + (u)[index+di*(1)+dj*(0)+dk*(1)]) + 8*((u)[index+di*(-1)+dj*(0)+dk*(2)] + (u)[index+di*(1)+dj*(0)+dk*(-2)] + (u)[index+di*(-2)+dj*(0)+dk*(1)] + (u)[index+di*(2)+dj*(0)+dk*(-1)]) - 8*((u)[index+di*(-1)+dj*(0)+dk*(-2)] + (u)[index+di*(1)+dj*(0)+dk*(2)] + (u)[index+di*(-2)+dj*(0)+dk*(-1)] + (u)[index+di*(2)+dj*(0)+dk*(1)]) + (u)[index+di*(-2)+dj*(0)+dk*(-2)] - (u)[index+di*(-2)+dj*(0)+dk*(2)] - (u)[index+di*(2)+dj*(0)+dk*(-2)] + (u)[index+di*(2)+dj*(0)+dk*(2)]))
-#define PDstandardNth21(u,i,j,k) (p1o144dxdy*(-64*((u)[index+di*(-1)+dj*(1)+dk*(0)] + (u)[index+di*(1)+dj*(-1)+dk*(0)]) + 64*((u)[index+di*(-1)+dj*(-1)+dk*(0)] + (u)[index+di*(1)+dj*(1)+dk*(0)]) + 8*((u)[index+di*(-1)+dj*(2)+dk*(0)] + (u)[index+di*(1)+dj*(-2)+dk*(0)] + (u)[index+di*(-2)+dj*(1)+dk*(0)] + (u)[index+di*(2)+dj*(-1)+dk*(0)]) - 8*((u)[index+di*(-1)+dj*(-2)+dk*(0)] + (u)[index+di*(1)+dj*(2)+dk*(0)] + (u)[index+di*(-2)+dj*(-1)+dk*(0)] + (u)[index+di*(2)+dj*(1)+dk*(0)]) + (u)[index+di*(-2)+dj*(-2)+dk*(0)] - (u)[index+di*(-2)+dj*(2)+dk*(0)] - (u)[index+di*(2)+dj*(-2)+dk*(0)] + (u)[index+di*(2)+dj*(2)+dk*(0)]))
-#define PDstandardNth23(u,i,j,k) (p1o144dydz*(-64*((u)[index+di*(0)+dj*(-1)+dk*(1)] + (u)[index+di*(0)+dj*(1)+dk*(-1)]) + 64*((u)[index+di*(0)+dj*(-1)+dk*(-1)] + (u)[index+di*(0)+dj*(1)+dk*(1)]) + 8*((u)[index+di*(0)+dj*(-1)+dk*(2)] + (u)[index+di*(0)+dj*(1)+dk*(-2)] + (u)[index+di*(0)+dj*(-2)+dk*(1)] + (u)[index+di*(0)+dj*(2)+dk*(-1)]) - 8*((u)[index+di*(0)+dj*(-1)+dk*(-2)] + (u)[index+di*(0)+dj*(1)+dk*(2)] + (u)[index+di*(0)+dj*(-2)+dk*(-1)] + (u)[index+di*(0)+dj*(2)+dk*(1)]) + (u)[index+di*(0)+dj*(-2)+dk*(-2)] - (u)[index+di*(0)+dj*(-2)+dk*(2)] - (u)[index+di*(0)+dj*(2)+dk*(-2)] + (u)[index+di*(0)+dj*(2)+dk*(2)]))
-#define PDstandardNth31(u,i,j,k) (p1o144dxdz*(-64*((u)[index+di*(-1)+dj*(0)+dk*(1)] + (u)[index+di*(1)+dj*(0)+dk*(-1)]) + 64*((u)[index+di*(-1)+dj*(0)+dk*(-1)] + (u)[index+di*(1)+dj*(0)+dk*(1)]) + 8*((u)[index+di*(-1)+dj*(0)+dk*(2)] + (u)[index+di*(1)+dj*(0)+dk*(-2)] + (u)[index+di*(-2)+dj*(0)+dk*(1)] + (u)[index+di*(2)+dj*(0)+dk*(-1)]) - 8*((u)[index+di*(-1)+dj*(0)+dk*(-2)] + (u)[index+di*(1)+dj*(0)+dk*(2)] + (u)[index+di*(-2)+dj*(0)+dk*(-1)] + (u)[index+di*(2)+dj*(0)+dk*(1)]) + (u)[index+di*(-2)+dj*(0)+dk*(-2)] - (u)[index+di*(-2)+dj*(0)+dk*(2)] - (u)[index+di*(2)+dj*(0)+dk*(-2)] + (u)[index+di*(2)+dj*(0)+dk*(2)]))
-#define PDstandardNth32(u,i,j,k) (p1o144dydz*(-64*((u)[index+di*(0)+dj*(-1)+dk*(1)] + (u)[index+di*(0)+dj*(1)+dk*(-1)]) + 64*((u)[index+di*(0)+dj*(-1)+dk*(-1)] + (u)[index+di*(0)+dj*(1)+dk*(1)]) + 8*((u)[index+di*(0)+dj*(-1)+dk*(2)] + (u)[index+di*(0)+dj*(1)+dk*(-2)] + (u)[index+di*(0)+dj*(-2)+dk*(1)] + (u)[index+di*(0)+dj*(2)+dk*(-1)]) - 8*((u)[index+di*(0)+dj*(-1)+dk*(-2)] + (u)[index+di*(0)+dj*(1)+dk*(2)] + (u)[index+di*(0)+dj*(-2)+dk*(-1)] + (u)[index+di*(0)+dj*(2)+dk*(1)]) + (u)[index+di*(0)+dj*(-2)+dk*(-2)] - (u)[index+di*(0)+dj*(-2)+dk*(2)] - (u)[index+di*(0)+dj*(2)+dk*(-2)] + (u)[index+di*(0)+dj*(2)+dk*(2)]))
+#define PDstandardNth1(u,i,j,k) (p1o12dx*(-8*vec_loadu((u)[index+di*(-1)+dj*(0)+dk*(0)]) + 8*vec_loadu((u)[index+di*(1)+dj*(0)+dk*(0)]) + vec_loadu((u)[index+di*(-2)+dj*(0)+dk*(0)]) - vec_loadu((u)[index+di*(2)+dj*(0)+dk*(0)])))
+#define PDstandardNth2(u,i,j,k) (p1o12dy*(-8*vec_loadu((u)[index+di*(0)+dj*(-1)+dk*(0)]) + 8*vec_loadu((u)[index+di*(0)+dj*(1)+dk*(0)]) + vec_loadu((u)[index+di*(0)+dj*(-2)+dk*(0)]) - vec_loadu((u)[index+di*(0)+dj*(2)+dk*(0)])))
+#define PDstandardNth3(u,i,j,k) (p1o12dz*(-8*vec_loadu((u)[index+di*(0)+dj*(0)+dk*(-1)]) + 8*vec_loadu((u)[index+di*(0)+dj*(0)+dk*(1)]) + vec_loadu((u)[index+di*(0)+dj*(0)+dk*(-2)]) - vec_loadu((u)[index+di*(0)+dj*(0)+dk*(2)])))
+#define PDstandardNth11(u,i,j,k) (pm1o12dx2*(30*vec_loadu((u)[index+di*(0)+dj*(0)+dk*(0)]) - 16*(vec_loadu((u)[index+di*(-1)+dj*(0)+dk*(0)]) + vec_loadu((u)[index+di*(1)+dj*(0)+dk*(0)])) + vec_loadu((u)[index+di*(-2)+dj*(0)+dk*(0)]) + vec_loadu((u)[index+di*(2)+dj*(0)+dk*(0)])))
+#define PDstandardNth22(u,i,j,k) (pm1o12dy2*(30*vec_loadu((u)[index+di*(0)+dj*(0)+dk*(0)]) - 16*(vec_loadu((u)[index+di*(0)+dj*(-1)+dk*(0)]) + vec_loadu((u)[index+di*(0)+dj*(1)+dk*(0)])) + vec_loadu((u)[index+di*(0)+dj*(-2)+dk*(0)]) + vec_loadu((u)[index+di*(0)+dj*(2)+dk*(0)])))
+#define PDstandardNth33(u,i,j,k) (pm1o12dz2*(30*vec_loadu((u)[index+di*(0)+dj*(0)+dk*(0)]) - 16*(vec_loadu((u)[index+di*(0)+dj*(0)+dk*(-1)]) + vec_loadu((u)[index+di*(0)+dj*(0)+dk*(1)])) + vec_loadu((u)[index+di*(0)+dj*(0)+dk*(-2)]) + vec_loadu((u)[index+di*(0)+dj*(0)+dk*(2)])))
+#define PDstandardNth12(u,i,j,k) (p1o144dxdy*(-64*(vec_loadu((u)[index+di*(-1)+dj*(1)+dk*(0)]) + vec_loadu((u)[index+di*(1)+dj*(-1)+dk*(0)])) + 64*(vec_loadu((u)[index+di*(-1)+dj*(-1)+dk*(0)]) + vec_loadu((u)[index+di*(1)+dj*(1)+dk*(0)])) + 8*(vec_loadu((u)[index+di*(-1)+dj*(2)+dk*(0)]) + vec_loadu((u)[index+di*(1)+dj*(-2)+dk*(0)]) + vec_loadu((u)[index+di*(-2)+dj*(1)+dk*(0)]) + vec_loadu((u)[index+di*(2)+dj*(-1)+dk*(0)])) - 8*(vec_loadu((u)[index+di*(-1)+dj*(-2)+dk*(0)]) + vec_loadu((u)[index+di*(1)+dj*(2)+dk*(0)]) + vec_loadu((u)[index+di*(-2)+dj*(-1)+dk*(0)]) + vec_loadu((u)[index+di*(2)+dj*(1)+dk*(0)])) + vec_loadu((u)[index+di*(-2)+dj*(-2)+dk*(0)]) - vec_loadu((u)[index+di*(-2)+dj*(2)+dk*(0)]) - vec_loadu((u)[index+di*(2)+dj*(-2)+dk*(0)]) + vec_loadu((u)[index+di*(2)+dj*(2)+dk*(0)])))
+#define PDstandardNth13(u,i,j,k) (p1o144dxdz*(-64*(vec_loadu((u)[index+di*(-1)+dj*(0)+dk*(1)]) + vec_loadu((u)[index+di*(1)+dj*(0)+dk*(-1)])) + 64*(vec_loadu((u)[index+di*(-1)+dj*(0)+dk*(-1)]) + vec_loadu((u)[index+di*(1)+dj*(0)+dk*(1)])) + 8*(vec_loadu((u)[index+di*(-1)+dj*(0)+dk*(2)]) + vec_loadu((u)[index+di*(1)+dj*(0)+dk*(-2)]) + vec_loadu((u)[index+di*(-2)+dj*(0)+dk*(1)]) + vec_loadu((u)[index+di*(2)+dj*(0)+dk*(-1)])) - 8*(vec_loadu((u)[index+di*(-1)+dj*(0)+dk*(-2)]) + vec_loadu((u)[index+di*(1)+dj*(0)+dk*(2)]) + vec_loadu((u)[index+di*(-2)+dj*(0)+dk*(-1)]) + vec_loadu((u)[index+di*(2)+dj*(0)+dk*(1)])) + vec_loadu((u)[index+di*(-2)+dj*(0)+dk*(-2)]) - vec_loadu((u)[index+di*(-2)+dj*(0)+dk*(2)]) - vec_loadu((u)[index+di*(2)+dj*(0)+dk*(-2)]) + vec_loadu((u)[index+di*(2)+dj*(0)+dk*(2)])))
+#define PDstandardNth21(u,i,j,k) (p1o144dxdy*(-64*(vec_loadu((u)[index+di*(-1)+dj*(1)+dk*(0)]) + vec_loadu((u)[index+di*(1)+dj*(-1)+dk*(0)])) + 64*(vec_loadu((u)[index+di*(-1)+dj*(-1)+dk*(0)]) + vec_loadu((u)[index+di*(1)+dj*(1)+dk*(0)])) + 8*(vec_loadu((u)[index+di*(-1)+dj*(2)+dk*(0)]) + vec_loadu((u)[index+di*(1)+dj*(-2)+dk*(0)]) + vec_loadu((u)[index+di*(-2)+dj*(1)+dk*(0)]) + vec_loadu((u)[index+di*(2)+dj*(-1)+dk*(0)])) - 8*(vec_loadu((u)[index+di*(-1)+dj*(-2)+dk*(0)]) + vec_loadu((u)[index+di*(1)+dj*(2)+dk*(0)]) + vec_loadu((u)[index+di*(-2)+dj*(-1)+dk*(0)]) + vec_loadu((u)[index+di*(2)+dj*(1)+dk*(0)])) + vec_loadu((u)[index+di*(-2)+dj*(-2)+dk*(0)]) - vec_loadu((u)[index+di*(-2)+dj*(2)+dk*(0)]) - vec_loadu((u)[index+di*(2)+dj*(-2)+dk*(0)]) + vec_loadu((u)[index+di*(2)+dj*(2)+dk*(0)])))
+#define PDstandardNth23(u,i,j,k) (p1o144dydz*(-64*(vec_loadu((u)[index+di*(0)+dj*(-1)+dk*(1)]) + vec_loadu((u)[index+di*(0)+dj*(1)+dk*(-1)])) + 64*(vec_loadu((u)[index+di*(0)+dj*(-1)+dk*(-1)]) + vec_loadu((u)[index+di*(0)+dj*(1)+dk*(1)])) + 8*(vec_loadu((u)[index+di*(0)+dj*(-1)+dk*(2)]) + vec_loadu((u)[index+di*(0)+dj*(1)+dk*(-2)]) + vec_loadu((u)[index+di*(0)+dj*(-2)+dk*(1)]) + vec_loadu((u)[index+di*(0)+dj*(2)+dk*(-1)])) - 8*(vec_loadu((u)[index+di*(0)+dj*(-1)+dk*(-2)]) + vec_loadu((u)[index+di*(0)+dj*(1)+dk*(2)]) + vec_loadu((u)[index+di*(0)+dj*(-2)+dk*(-1)]) + vec_loadu((u)[index+di*(0)+dj*(2)+dk*(1)])) + vec_loadu((u)[index+di*(0)+dj*(-2)+dk*(-2)]) - vec_loadu((u)[index+di*(0)+dj*(-2)+dk*(2)]) - vec_loadu((u)[index+di*(0)+dj*(2)+dk*(-2)]) + vec_loadu((u)[index+di*(0)+dj*(2)+dk*(2)])))
+#define PDstandardNth31(u,i,j,k) (p1o144dxdz*(-64*(vec_loadu((u)[index+di*(-1)+dj*(0)+dk*(1)]) + vec_loadu((u)[index+di*(1)+dj*(0)+dk*(-1)])) + 64*(vec_loadu((u)[index+di*(-1)+dj*(0)+dk*(-1)]) + vec_loadu((u)[index+di*(1)+dj*(0)+dk*(1)])) + 8*(vec_loadu((u)[index+di*(-1)+dj*(0)+dk*(2)]) + vec_loadu((u)[index+di*(1)+dj*(0)+dk*(-2)]) + vec_loadu((u)[index+di*(-2)+dj*(0)+dk*(1)]) + vec_loadu((u)[index+di*(2)+dj*(0)+dk*(-1)])) - 8*(vec_loadu((u)[index+di*(-1)+dj*(0)+dk*(-2)]) + vec_loadu((u)[index+di*(1)+dj*(0)+dk*(2)]) + vec_loadu((u)[index+di*(-2)+dj*(0)+dk*(-1)]) + vec_loadu((u)[index+di*(2)+dj*(0)+dk*(1)])) + vec_loadu((u)[index+di*(-2)+dj*(0)+dk*(-2)]) - vec_loadu((u)[index+di*(-2)+dj*(0)+dk*(2)]) - vec_loadu((u)[index+di*(2)+dj*(0)+dk*(-2)]) + vec_loadu((u)[index+di*(2)+dj*(0)+dk*(2)])))
+#define PDstandardNth32(u,i,j,k) (p1o144dydz*(-64*(vec_loadu((u)[index+di*(0)+dj*(-1)+dk*(1)]) + vec_loadu((u)[index+di*(0)+dj*(1)+dk*(-1)])) + 64*(vec_loadu((u)[index+di*(0)+dj*(-1)+dk*(-1)]) + vec_loadu((u)[index+di*(0)+dj*(1)+dk*(1)])) + 8*(vec_loadu((u)[index+di*(0)+dj*(-1)+dk*(2)]) + vec_loadu((u)[index+di*(0)+dj*(1)+dk*(-2)]) + vec_loadu((u)[index+di*(0)+dj*(-2)+dk*(1)]) + vec_loadu((u)[index+di*(0)+dj*(2)+dk*(-1)])) - 8*(vec_loadu((u)[index+di*(0)+dj*(-1)+dk*(-2)]) + vec_loadu((u)[index+di*(0)+dj*(1)+dk*(2)]) + vec_loadu((u)[index+di*(0)+dj*(-2)+dk*(-1)]) + vec_loadu((u)[index+di*(0)+dj*(2)+dk*(1)])) + vec_loadu((u)[index+di*(0)+dj*(-2)+dk*(-2)]) - vec_loadu((u)[index+di*(0)+dj*(-2)+dk*(2)]) - vec_loadu((u)[index+di*(0)+dj*(2)+dk*(-2)]) + vec_loadu((u)[index+di*(0)+dj*(2)+dk*(2)])))
diff --git a/ML_WaveToyFO/src/WTFO_Gaussian.c b/ML_WaveToyFO/src/WTFO_Gaussian.cc
index 85a4f97..dd7821a 100644
--- a/ML_WaveToyFO/src/WTFO_Gaussian.c
+++ b/ML_WaveToyFO/src/WTFO_Gaussian.cc
@@ -11,6 +11,7 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "Vectors.hh"
 #include "loopcontrol.h"
 
 /* Define macros used in calculations */
@@ -20,7 +21,7 @@
 #define CUB(x) ((x) * (x) * (x))
 #define QAD(x) ((x) * (x) * (x) * (x))
 
-void WTFO_Gaussian_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void WTFO_Gaussian_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
@@ -28,17 +29,6 @@ void WTFO_Gaussian_Body(cGH const * restrict const cctkGH, int const dir, int co
   
   /* Declare finite differencing variables */
   
-  /* Declare predefined quantities */
-  // CCTK_REAL p1o12dx = INITVALUE;
-  // CCTK_REAL p1o12dy = INITVALUE;
-  // CCTK_REAL p1o12dz = INITVALUE;
-  // CCTK_REAL p1o144dxdy = INITVALUE;
-  // CCTK_REAL p1o144dxdz = INITVALUE;
-  // CCTK_REAL p1o144dydz = INITVALUE;
-  // CCTK_REAL pm1o12dx2 = INITVALUE;
-  // CCTK_REAL pm1o12dy2 = INITVALUE;
-  // CCTK_REAL pm1o12dz2 = INITVALUE;
-  
   if (verbose > 1)
   {
     CCTK_VInfo(CCTK_THORNSTRING,"Entering WTFO_Gaussian_Body");
@@ -88,55 +78,39 @@ void WTFO_Gaussian_Body(cGH const * restrict const cctkGH, int const dir, int co
             cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
     // int index = INITVALUE;
-    // int subblock_index = INITVALUE;
     int const index = CCTK_GFINDEX3D(cctkGH,i,j,k);
-    int const subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2]));
-    
-    /* Declare shorthands */
-    
-    /* Declare local copies of grid functions */
-    // CCTK_REAL rhoL = INITVALUE;
-    // CCTK_REAL uL = INITVALUE;
-    // CCTK_REAL v1L = INITVALUE, v2L = INITVALUE, v3L = INITVALUE;
-    /* Declare precomputed derivatives*/
-    
     /* Declare derivatives */
     
     /* Assign local copies of grid functions */
     
-    /* Assign local copies of subblock grid functions */
-    
     /* Include user supplied include files */
     
-    /* Precompute derivatives (new style) */
-    
-    /* Precompute derivatives (old style) */
+    /* Precompute derivatives */
     
     /* Calculate temporaries and grid functions */
-    CCTK_REAL const uL  =  0;
+    CCTK_REAL_VEC uL = 0;
     
-    CCTK_REAL const v1L  =  0;
+    CCTK_REAL_VEC v1L = 0;
     
-    CCTK_REAL const v2L  =  0;
+    CCTK_REAL_VEC v2L = 0;
     
-    CCTK_REAL const v3L  =  0;
+    CCTK_REAL_VEC v3L = 0;
     
-    CCTK_REAL const rhoL  =  0;
+    CCTK_REAL_VEC rhoL = 0;
     
     
     /* Copy local copies back to grid functions */
-    rho[index] = rhoL;
-    u[index] = uL;
-    v1[index] = v1L;
-    v2[index] = v2L;
-    v3[index] = v3L;
-    
-    /* Copy local copies back to subblock grid functions */
+    vec_store_nta(rho[index],rhoL);
+    vec_store_nta(u[index],uL);
+    vec_store_nta(v1[index],v1L);
+    vec_store_nta(v2[index],v2L);
+    vec_store_nta(v3[index],v3L);
+  i += CCTK_REAL_VEC_SIZE-1;
   }
   LC_ENDLOOP3 (WTFO_Gaussian);
 }
 
-void WTFO_Gaussian(CCTK_ARGUMENTS)
+extern "C" void WTFO_Gaussian(CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
diff --git a/ML_WaveToyFO/src/WTFO_RHS.c b/ML_WaveToyFO/src/WTFO_RHS.cc
index 6c0477c..9961c5b 100644
--- a/ML_WaveToyFO/src/WTFO_RHS.c
+++ b/ML_WaveToyFO/src/WTFO_RHS.cc
@@ -11,6 +11,7 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "Vectors.hh"
 #include "loopcontrol.h"
 
 /* Define macros used in calculations */
@@ -20,7 +21,7 @@
 #define CUB(x) ((x) * (x) * (x))
 #define QAD(x) ((x) * (x) * (x) * (x))
 
-void WTFO_RHS_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void WTFO_RHS_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
@@ -28,17 +29,6 @@ void WTFO_RHS_Body(cGH const * restrict const cctkGH, int const dir, int const f
   
   /* Declare finite differencing variables */
   
-  /* Declare predefined quantities */
-  // CCTK_REAL p1o12dx = INITVALUE;
-  // CCTK_REAL p1o12dy = INITVALUE;
-  // CCTK_REAL p1o12dz = INITVALUE;
-  // CCTK_REAL p1o144dxdy = INITVALUE;
-  // CCTK_REAL p1o144dxdz = INITVALUE;
-  // CCTK_REAL p1o144dydz = INITVALUE;
-  // CCTK_REAL pm1o12dx2 = INITVALUE;
-  // CCTK_REAL pm1o12dy2 = INITVALUE;
-  // CCTK_REAL pm1o12dz2 = INITVALUE;
-  
   if (verbose > 1)
   {
     CCTK_VInfo(CCTK_THORNSTRING,"Entering WTFO_RHS_Body");
@@ -88,71 +78,56 @@ void WTFO_RHS_Body(cGH const * restrict const cctkGH, int const dir, int const f
             cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
     // int index = INITVALUE;
-    // int subblock_index = INITVALUE;
     int const index = CCTK_GFINDEX3D(cctkGH,i,j,k);
-    int const subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2]));
-    
-    /* Declare shorthands */
-    
-    /* Declare local copies of grid functions */
-    // CCTK_REAL rhoL = INITVALUE, rhorhsL = INITVALUE;
-    // CCTK_REAL urhsL = INITVALUE;
-    // CCTK_REAL v1L = INITVALUE, v1rhsL = INITVALUE, v2L = INITVALUE, v2rhsL = INITVALUE, v3L = INITVALUE, v3rhsL = INITVALUE;
-    /* Declare precomputed derivatives*/
-    
     /* Declare derivatives */
-    // CCTK_REAL PDstandardNth1rho = INITVALUE;
-    // CCTK_REAL PDstandardNth2rho = INITVALUE;
-    // CCTK_REAL PDstandardNth3rho = INITVALUE;
-    // CCTK_REAL PDstandardNth1v1 = INITVALUE;
-    // CCTK_REAL PDstandardNth2v2 = INITVALUE;
-    // CCTK_REAL PDstandardNth3v3 = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth1rho = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth2rho = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth3rho = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth1v1 = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth2v2 = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth3v3 = INITVALUE;
     
     /* Assign local copies of grid functions */
-    CCTK_REAL const rhoL = rho[index];
-    CCTK_REAL const v1L = v1[index];
-    CCTK_REAL const v2L = v2[index];
-    CCTK_REAL const v3L = v3[index];
-    
-    /* Assign local copies of subblock grid functions */
+    CCTK_REAL_VEC  rhoL = vec_load(rho[index]);
+    CCTK_REAL_VEC  v1L = vec_load(v1[index]);
+    CCTK_REAL_VEC  v2L = vec_load(v2[index]);
+    CCTK_REAL_VEC  v3L = vec_load(v3[index]);
     
     /* Include user supplied include files */
     
-    /* Precompute derivatives (new style) */
-    CCTK_REAL const PDstandardNth1rho = PDstandardNth1(rho, i, j, k);
-    CCTK_REAL const PDstandardNth2rho = PDstandardNth2(rho, i, j, k);
-    CCTK_REAL const PDstandardNth3rho = PDstandardNth3(rho, i, j, k);
-    CCTK_REAL const PDstandardNth1v1 = PDstandardNth1(v1, i, j, k);
-    CCTK_REAL const PDstandardNth2v2 = PDstandardNth2(v2, i, j, k);
-    CCTK_REAL const PDstandardNth3v3 = PDstandardNth3(v3, i, j, k);
-    
-    /* Precompute derivatives (old style) */
+    /* Precompute derivatives */
+    CCTK_REAL_VEC const PDstandardNth1rho = PDstandardNth1(rho, i, j, k);
+    CCTK_REAL_VEC const PDstandardNth2rho = PDstandardNth2(rho, i, j, k);
+    CCTK_REAL_VEC const PDstandardNth3rho = PDstandardNth3(rho, i, j, k);
+    CCTK_REAL_VEC const PDstandardNth1v1 = PDstandardNth1(v1, i, j, k);
+    CCTK_REAL_VEC const PDstandardNth2v2 = PDstandardNth2(v2, i, j, k);
+    CCTK_REAL_VEC const PDstandardNth3v3 = PDstandardNth3(v3, i, j, k);
     
     /* Calculate temporaries and grid functions */
-    CCTK_REAL const urhsL  =  rhoL;
+    CCTK_REAL_VEC urhsL = rhoL;
     
-    CCTK_REAL const rhorhsL  =  PDstandardNth1v1 + PDstandardNth2v2 + PDstandardNth3v3;
+    CCTK_REAL_VEC rhorhsL = PDstandardNth1v1 + PDstandardNth2v2 + 
+      PDstandardNth3v3;
     
-    CCTK_REAL const v1rhsL  =  PDstandardNth1rho;
+    CCTK_REAL_VEC v1rhsL = PDstandardNth1rho;
     
-    CCTK_REAL const v2rhsL  =  PDstandardNth2rho;
+    CCTK_REAL_VEC v2rhsL = PDstandardNth2rho;
     
-    CCTK_REAL const v3rhsL  =  PDstandardNth3rho;
+    CCTK_REAL_VEC v3rhsL = PDstandardNth3rho;
     
     
     /* Copy local copies back to grid functions */
-    rhorhs[index] = rhorhsL;
-    urhs[index] = urhsL;
-    v1rhs[index] = v1rhsL;
-    v2rhs[index] = v2rhsL;
-    v3rhs[index] = v3rhsL;
-    
-    /* Copy local copies back to subblock grid functions */
+    vec_store_nta(rhorhs[index],rhorhsL);
+    vec_store_nta(urhs[index],urhsL);
+    vec_store_nta(v1rhs[index],v1rhsL);
+    vec_store_nta(v2rhs[index],v2rhsL);
+    vec_store_nta(v3rhs[index],v3rhsL);
+  i += CCTK_REAL_VEC_SIZE-1;
   }
   LC_ENDLOOP3 (WTFO_RHS);
 }
 
-void WTFO_RHS(CCTK_ARGUMENTS)
+extern "C" void WTFO_RHS(CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
diff --git a/ML_WaveToyFO/src/WTFO_constraints.c b/ML_WaveToyFO/src/WTFO_constraints.cc
index b865ffc..bdb4471 100644
--- a/ML_WaveToyFO/src/WTFO_constraints.c
+++ b/ML_WaveToyFO/src/WTFO_constraints.cc
@@ -11,6 +11,7 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "Vectors.hh"
 #include "loopcontrol.h"
 
 /* Define macros used in calculations */
@@ -20,7 +21,7 @@
 #define CUB(x) ((x) * (x) * (x))
 #define QAD(x) ((x) * (x) * (x) * (x))
 
-void WTFO_constraints_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void WTFO_constraints_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
@@ -28,17 +29,6 @@ void WTFO_constraints_Body(cGH const * restrict const cctkGH, int const dir, int
   
   /* Declare finite differencing variables */
   
-  /* Declare predefined quantities */
-  // CCTK_REAL p1o12dx = INITVALUE;
-  // CCTK_REAL p1o12dy = INITVALUE;
-  // CCTK_REAL p1o12dz = INITVALUE;
-  // CCTK_REAL p1o144dxdy = INITVALUE;
-  // CCTK_REAL p1o144dxdz = INITVALUE;
-  // CCTK_REAL p1o144dydz = INITVALUE;
-  // CCTK_REAL pm1o12dx2 = INITVALUE;
-  // CCTK_REAL pm1o12dy2 = INITVALUE;
-  // CCTK_REAL pm1o12dz2 = INITVALUE;
-  
   if (verbose > 1)
   {
     CCTK_VInfo(CCTK_THORNSTRING,"Entering WTFO_constraints_Body");
@@ -88,63 +78,48 @@ void WTFO_constraints_Body(cGH const * restrict const cctkGH, int const dir, int
             cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
     // int index = INITVALUE;
-    // int subblock_index = INITVALUE;
     int const index = CCTK_GFINDEX3D(cctkGH,i,j,k);
-    int const subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2]));
-    
-    /* Declare shorthands */
-    
-    /* Declare local copies of grid functions */
-    // CCTK_REAL v1L = INITVALUE, v2L = INITVALUE, v3L = INITVALUE;
-    // CCTK_REAL w1L = INITVALUE, w2L = INITVALUE, w3L = INITVALUE;
-    /* Declare precomputed derivatives*/
-    
     /* Declare derivatives */
-    // CCTK_REAL PDstandardNth2v1 = INITVALUE;
-    // CCTK_REAL PDstandardNth3v1 = INITVALUE;
-    // CCTK_REAL PDstandardNth1v2 = INITVALUE;
-    // CCTK_REAL PDstandardNth3v2 = INITVALUE;
-    // CCTK_REAL PDstandardNth1v3 = INITVALUE;
-    // CCTK_REAL PDstandardNth2v3 = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth2v1 = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth3v1 = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth1v2 = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth3v2 = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth1v3 = INITVALUE;
+    // CCTK_REAL_VEC PDstandardNth2v3 = INITVALUE;
     
     /* Assign local copies of grid functions */
-    CCTK_REAL const v1L = v1[index];
-    CCTK_REAL const v2L = v2[index];
-    CCTK_REAL const v3L = v3[index];
-    
-    /* Assign local copies of subblock grid functions */
+    CCTK_REAL_VEC  v1L = vec_load(v1[index]);
+    CCTK_REAL_VEC  v2L = vec_load(v2[index]);
+    CCTK_REAL_VEC  v3L = vec_load(v3[index]);
     
     /* Include user supplied include files */
     
-    /* Precompute derivatives (new style) */
-    CCTK_REAL const PDstandardNth2v1 = PDstandardNth2(v1, i, j, k);
-    CCTK_REAL const PDstandardNth3v1 = PDstandardNth3(v1, i, j, k);
-    CCTK_REAL const PDstandardNth1v2 = PDstandardNth1(v2, i, j, k);
-    CCTK_REAL const PDstandardNth3v2 = PDstandardNth3(v2, i, j, k);
-    CCTK_REAL const PDstandardNth1v3 = PDstandardNth1(v3, i, j, k);
-    CCTK_REAL const PDstandardNth2v3 = PDstandardNth2(v3, i, j, k);
-    
-    /* Precompute derivatives (old style) */
+    /* Precompute derivatives */
+    CCTK_REAL_VEC const PDstandardNth2v1 = PDstandardNth2(v1, i, j, k);
+    CCTK_REAL_VEC const PDstandardNth3v1 = PDstandardNth3(v1, i, j, k);
+    CCTK_REAL_VEC const PDstandardNth1v2 = PDstandardNth1(v2, i, j, k);
+    CCTK_REAL_VEC const PDstandardNth3v2 = PDstandardNth3(v2, i, j, k);
+    CCTK_REAL_VEC const PDstandardNth1v3 = PDstandardNth1(v3, i, j, k);
+    CCTK_REAL_VEC const PDstandardNth2v3 = PDstandardNth2(v3, i, j, k);
     
     /* Calculate temporaries and grid functions */
-    CCTK_REAL const w1L  =  -PDstandardNth2v3 + PDstandardNth3v2;
+    CCTK_REAL_VEC w1L = -PDstandardNth2v3 + PDstandardNth3v2;
     
-    CCTK_REAL const w2L  =  PDstandardNth1v3 - PDstandardNth3v1;
+    CCTK_REAL_VEC w2L = PDstandardNth1v3 - PDstandardNth3v1;
     
-    CCTK_REAL const w3L  =  -PDstandardNth1v2 + PDstandardNth2v1;
+    CCTK_REAL_VEC w3L = -PDstandardNth1v2 + PDstandardNth2v1;
     
     
     /* Copy local copies back to grid functions */
-    w1[index] = w1L;
-    w2[index] = w2L;
-    w3[index] = w3L;
-    
-    /* Copy local copies back to subblock grid functions */
+    vec_store_nta(w1[index],w1L);
+    vec_store_nta(w2[index],w2L);
+    vec_store_nta(w3[index],w3L);
+  i += CCTK_REAL_VEC_SIZE-1;
   }
   LC_ENDLOOP3 (WTFO_constraints);
 }
 
-void WTFO_constraints(CCTK_ARGUMENTS)
+extern "C" void WTFO_constraints(CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
diff --git a/ML_WaveToyFO/src/make.code.defn b/ML_WaveToyFO/src/make.code.defn
index a18c783..a51a3f3 100644
--- a/ML_WaveToyFO/src/make.code.defn
+++ b/ML_WaveToyFO/src/make.code.defn
@@ -1,3 +1,3 @@
 # File produced by Kranc
 
-SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c WTFO_Gaussian.c WTFO_RHS.c WTFO_constraints.c Boundaries.c 
+SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c WTFO_Gaussian.cc WTFO_RHS.cc WTFO_constraints.cc Boundaries.c