Regenerate examples (except Advect)

Advect doesn't work at the moment due to a problem with TensorTools.

10 files changed, 806 insertions, 280 deletions

diff --git a/Examples/EM/configuration.ccl b/Examples/EM/configuration.ccl
index 023aac7..0a66ec2 100644
--- a/Examples/EM/configuration.ccl
+++ b/Examples/EM/configuration.ccl
@@ -1,3 +1,6 @@
 # File produced by Kranc
 
 REQUIRES GenericFD
+OPTIONAL LoopControl
+{
+}
diff --git a/Examples/EM/interface.ccl b/Examples/EM/interface.ccl
index 58b38f2..46d63e1 100644
--- a/Examples/EM/interface.ccl
+++ b/Examples/EM/interface.ccl
@@ -10,6 +10,7 @@ USES INCLUDE: GenericFD.h
 USES INCLUDE: Symmetry.h
 USES INCLUDE: sbp_calc_coeffs.h
 USES INCLUDE: Boundary.h
+USES INCLUDE: loopcontrol.h
 
 CCTK_INT FUNCTION MoLRegisterEvolved(CCTK_INT IN EvolvedIndex, CCTK_INT IN RHSIndex)
 USES FUNCTION MoLRegisterEvolved
diff --git a/Examples/EM/param.ccl b/Examples/EM/param.ccl
index fad7962..4d95cbd 100644
--- a/Examples/EM/param.ccl
+++ b/Examples/EM/param.ccl
@@ -8,6 +8,7 @@ shares: GenericFD
 shares: MethodOfLines
 
 USES CCTK_INT MoL_Num_Evolved_Vars
+USES CCTK_INT MoL_Num_ArrayEvolved_Vars
 
 restricted:
 CCTK_INT verbose "verbose" STEERABLE=ALWAYS
@@ -34,6 +35,12 @@ CCTK_INT EM_MaxNumEvolvedVars "Number of evolved variables used by this thorn" A
 } 6
 
 restricted:
+CCTK_INT EM_MaxNumArrayEvolvedVars "Number of Array evolved variables used by this thorn" ACCUMULATOR-BASE=MethodofLines::MoL_Num_ArrayEvolved_Vars STEERABLE=RECOVER
+{
+  0:0 :: "Number of Array evolved variables used by this thorn"
+} 0
+
+restricted:
 CCTK_INT timelevels "Number of active timelevels" STEERABLE=RECOVER
 {
   0:3 :: ""
@@ -46,6 +53,12 @@ CCTK_INT rhs_timelevels "Number of active RHS timelevels" STEERABLE=RECOVER
 } 1
 
 restricted:
+CCTK_INT other_timelevels "Number of active timelevels for non-evolved grid functions" STEERABLE=RECOVER
+{
+  0:3 :: ""
+} 1
+
+restricted:
 CCTK_INT EM_initial_calc_every "EM_initial_calc_every" STEERABLE=ALWAYS
 {
   *:* :: ""
diff --git a/Examples/EM/schedule.ccl b/Examples/EM/schedule.ccl
index fea3ede..3dfb80d 100644
--- a/Examples/EM/schedule.ccl
+++ b/Examples/EM/schedule.ccl
@@ -1,9 +1,15 @@
 # File produced by Kranc
 
 
-STORAGE: constraints[1]
+if (other_timelevels == 1)
+{
+  STORAGE: constraints[1]
+}
 
-STORAGE: endens[1]
+if (other_timelevels == 1)
+{
+  STORAGE: endens[1]
+}
 
 if (timelevels == 1)
 {
@@ -95,6 +101,28 @@ schedule EM_constraints in EM_constraints_group
   LANG: C
 } "EM_constraints"
 
+schedule EM_constraints_SelectBCs in EM_constraints_bc_group
+{
+  LANG: C
+  OPTIONS: level
+  SYNC: constraints
+} "EM_constraints_SelectBCs"
+
+schedule group ApplyBCs as EM_constraints_ApplyBCs in EM_constraints_bc_group after EM_constraints_SelectBCs
+{
+  # no language specified
+} "Apply BCs for groups set in EM_constraints"
+
+schedule group EM_constraints_bc_group in EM_constraints_group after EM_constraints
+{
+  # no language specified
+} "EM_constraints"
+
+schedule group EM_constraints_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
+{
+  # no language specified
+} "EM_constraints"
+
 schedule group EM_energy_group in MoL_PseudoEvolution
 {
   # no language specified
diff --git a/Examples/EM/src/Differencing.h b/Examples/EM/src/Differencing.h
index 7740e61..e6b9ac7 100644
--- a/Examples/EM/src/Differencing.h
+++ b/Examples/EM/src/Differencing.h
@@ -1,36 +1,432 @@
-#define PDstandard2nd1(u) (p1o2dx*(-(u)[di*(-1)+dj*(0)+dk*(0)] + (u)[di*(1)+dj*(0)+dk*(0)]))
-#define PDstandard2nd2(u) (p1o2dy*(-(u)[di*(0)+dj*(-1)+dk*(0)] + (u)[di*(0)+dj*(1)+dk*(0)]))
-#define PDstandard2nd3(u) (p1o2dz*(-(u)[di*(0)+dj*(0)+dk*(-1)] + (u)[di*(0)+dj*(0)+dk*(1)]))
-#define PDstandard2nd11(u) (p1odx2*(-2*(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(-1)+dj*(0)+dk*(0)] + (u)[di*(1)+dj*(0)+dk*(0)]))
-#define PDstandard2nd22(u) (p1ody2*(-2*(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(0)+dj*(-1)+dk*(0)] + (u)[di*(0)+dj*(1)+dk*(0)]))
-#define PDstandard2nd33(u) (p1odz2*(-2*(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(0)+dj*(0)+dk*(-1)] + (u)[di*(0)+dj*(0)+dk*(1)]))
-#define PDstandard2nd12(u) (p1o4dxdy*((u)[di*(-1)+dj*(-1)+dk*(0)] - (u)[di*(-1)+dj*(1)+dk*(0)] - (u)[di*(1)+dj*(-1)+dk*(0)] + (u)[di*(1)+dj*(1)+dk*(0)]))
-#define PDstandard2nd13(u) (p1o4dxdz*((u)[di*(-1)+dj*(0)+dk*(-1)] - (u)[di*(-1)+dj*(0)+dk*(1)] - (u)[di*(1)+dj*(0)+dk*(-1)] + (u)[di*(1)+dj*(0)+dk*(1)]))
-#define PDstandard2nd21(u) (p1o4dxdy*((u)[di*(-1)+dj*(-1)+dk*(0)] - (u)[di*(-1)+dj*(1)+dk*(0)] - (u)[di*(1)+dj*(-1)+dk*(0)] + (u)[di*(1)+dj*(1)+dk*(0)]))
-#define PDstandard2nd23(u) (p1o4dydz*((u)[di*(0)+dj*(-1)+dk*(-1)] - (u)[di*(0)+dj*(-1)+dk*(1)] - (u)[di*(0)+dj*(1)+dk*(-1)] + (u)[di*(0)+dj*(1)+dk*(1)]))
-#define PDstandard2nd31(u) (p1o4dxdz*((u)[di*(-1)+dj*(0)+dk*(-1)] - (u)[di*(-1)+dj*(0)+dk*(1)] - (u)[di*(1)+dj*(0)+dk*(-1)] + (u)[di*(1)+dj*(0)+dk*(1)]))
-#define PDstandard2nd32(u) (p1o4dydz*((u)[di*(0)+dj*(-1)+dk*(-1)] - (u)[di*(0)+dj*(-1)+dk*(1)] - (u)[di*(0)+dj*(1)+dk*(-1)] + (u)[di*(0)+dj*(1)+dk*(1)]))
-#define PDstandard4th1(u) (p1o12dx*(-8*(u)[di*(-1)+dj*(0)+dk*(0)] + 8*(u)[di*(1)+dj*(0)+dk*(0)] + (u)[di*(-2)+dj*(0)+dk*(0)] - (u)[di*(2)+dj*(0)+dk*(0)]))
-#define PDstandard4th2(u) (p1o12dy*(-8*(u)[di*(0)+dj*(-1)+dk*(0)] + 8*(u)[di*(0)+dj*(1)+dk*(0)] + (u)[di*(0)+dj*(-2)+dk*(0)] - (u)[di*(0)+dj*(2)+dk*(0)]))
-#define PDstandard4th3(u) (p1o12dz*(-8*(u)[di*(0)+dj*(0)+dk*(-1)] + 8*(u)[di*(0)+dj*(0)+dk*(1)] + (u)[di*(0)+dj*(0)+dk*(-2)] - (u)[di*(0)+dj*(0)+dk*(2)]))
-#define PDstandard4th11(u) (pm1o12dx2*(30*(u)[di*(0)+dj*(0)+dk*(0)] - 16*((u)[di*(-1)+dj*(0)+dk*(0)] + (u)[di*(1)+dj*(0)+dk*(0)]) + (u)[di*(-2)+dj*(0)+dk*(0)] + (u)[di*(2)+dj*(0)+dk*(0)]))
-#define PDstandard4th22(u) (pm1o12dy2*(30*(u)[di*(0)+dj*(0)+dk*(0)] - 16*((u)[di*(0)+dj*(-1)+dk*(0)] + (u)[di*(0)+dj*(1)+dk*(0)]) + (u)[di*(0)+dj*(-2)+dk*(0)] + (u)[di*(0)+dj*(2)+dk*(0)]))
-#define PDstandard4th33(u) (pm1o12dz2*(30*(u)[di*(0)+dj*(0)+dk*(0)] - 16*((u)[di*(0)+dj*(0)+dk*(-1)] + (u)[di*(0)+dj*(0)+dk*(1)]) + (u)[di*(0)+dj*(0)+dk*(-2)] + (u)[di*(0)+dj*(0)+dk*(2)]))
-#define PDstandard4th12(u) (p1o144dxdy*(-64*((u)[di*(-1)+dj*(1)+dk*(0)] + (u)[di*(1)+dj*(-1)+dk*(0)]) + 64*((u)[di*(-1)+dj*(-1)+dk*(0)] + (u)[di*(1)+dj*(1)+dk*(0)]) + 8*((u)[di*(-1)+dj*(2)+dk*(0)] + (u)[di*(1)+dj*(-2)+dk*(0)] + (u)[di*(-2)+dj*(1)+dk*(0)] + (u)[di*(2)+dj*(-1)+dk*(0)]) - 8*((u)[di*(-1)+dj*(-2)+dk*(0)] + (u)[di*(1)+dj*(2)+dk*(0)] + (u)[di*(-2)+dj*(-1)+dk*(0)] + (u)[di*(2)+dj*(1)+dk*(0)]) + (u)[di*(-2)+dj*(-2)+dk*(0)] - (u)[di*(-2)+dj*(2)+dk*(0)] - (u)[di*(2)+dj*(-2)+dk*(0)] + (u)[di*(2)+dj*(2)+dk*(0)]))
-#define PDstandard4th13(u) (p1o144dxdz*(-64*((u)[di*(-1)+dj*(0)+dk*(1)] + (u)[di*(1)+dj*(0)+dk*(-1)]) + 64*((u)[di*(-1)+dj*(0)+dk*(-1)] + (u)[di*(1)+dj*(0)+dk*(1)]) + 8*((u)[di*(-1)+dj*(0)+dk*(2)] + (u)[di*(1)+dj*(0)+dk*(-2)] + (u)[di*(-2)+dj*(0)+dk*(1)] + (u)[di*(2)+dj*(0)+dk*(-1)]) - 8*((u)[di*(-1)+dj*(0)+dk*(-2)] + (u)[di*(1)+dj*(0)+dk*(2)] + (u)[di*(-2)+dj*(0)+dk*(-1)] + (u)[di*(2)+dj*(0)+dk*(1)]) + (u)[di*(-2)+dj*(0)+dk*(-2)] - (u)[di*(-2)+dj*(0)+dk*(2)] - (u)[di*(2)+dj*(0)+dk*(-2)] + (u)[di*(2)+dj*(0)+dk*(2)]))
-#define PDstandard4th21(u) (p1o144dxdy*(-64*((u)[di*(-1)+dj*(1)+dk*(0)] + (u)[di*(1)+dj*(-1)+dk*(0)]) + 64*((u)[di*(-1)+dj*(-1)+dk*(0)] + (u)[di*(1)+dj*(1)+dk*(0)]) + 8*((u)[di*(-1)+dj*(2)+dk*(0)] + (u)[di*(1)+dj*(-2)+dk*(0)] + (u)[di*(-2)+dj*(1)+dk*(0)] + (u)[di*(2)+dj*(-1)+dk*(0)]) - 8*((u)[di*(-1)+dj*(-2)+dk*(0)] + (u)[di*(1)+dj*(2)+dk*(0)] + (u)[di*(-2)+dj*(-1)+dk*(0)] + (u)[di*(2)+dj*(1)+dk*(0)]) + (u)[di*(-2)+dj*(-2)+dk*(0)] - (u)[di*(-2)+dj*(2)+dk*(0)] - (u)[di*(2)+dj*(-2)+dk*(0)] + (u)[di*(2)+dj*(2)+dk*(0)]))
-#define PDstandard4th23(u) (p1o144dydz*(-64*((u)[di*(0)+dj*(-1)+dk*(1)] + (u)[di*(0)+dj*(1)+dk*(-1)]) + 64*((u)[di*(0)+dj*(-1)+dk*(-1)] + (u)[di*(0)+dj*(1)+dk*(1)]) + 8*((u)[di*(0)+dj*(-1)+dk*(2)] + (u)[di*(0)+dj*(1)+dk*(-2)] + (u)[di*(0)+dj*(-2)+dk*(1)] + (u)[di*(0)+dj*(2)+dk*(-1)]) - 8*((u)[di*(0)+dj*(-1)+dk*(-2)] + (u)[di*(0)+dj*(1)+dk*(2)] + (u)[di*(0)+dj*(-2)+dk*(-1)] + (u)[di*(0)+dj*(2)+dk*(1)]) + (u)[di*(0)+dj*(-2)+dk*(-2)] - (u)[di*(0)+dj*(-2)+dk*(2)] - (u)[di*(0)+dj*(2)+dk*(-2)] + (u)[di*(0)+dj*(2)+dk*(2)]))
-#define PDstandard4th31(u) (p1o144dxdz*(-64*((u)[di*(-1)+dj*(0)+dk*(1)] + (u)[di*(1)+dj*(0)+dk*(-1)]) + 64*((u)[di*(-1)+dj*(0)+dk*(-1)] + (u)[di*(1)+dj*(0)+dk*(1)]) + 8*((u)[di*(-1)+dj*(0)+dk*(2)] + (u)[di*(1)+dj*(0)+dk*(-2)] + (u)[di*(-2)+dj*(0)+dk*(1)] + (u)[di*(2)+dj*(0)+dk*(-1)]) - 8*((u)[di*(-1)+dj*(0)+dk*(-2)] + (u)[di*(1)+dj*(0)+dk*(2)] + (u)[di*(-2)+dj*(0)+dk*(-1)] + (u)[di*(2)+dj*(0)+dk*(1)]) + (u)[di*(-2)+dj*(0)+dk*(-2)] - (u)[di*(-2)+dj*(0)+dk*(2)] - (u)[di*(2)+dj*(0)+dk*(-2)] + (u)[di*(2)+dj*(0)+dk*(2)]))
-#define PDstandard4th32(u) (p1o144dydz*(-64*((u)[di*(0)+dj*(-1)+dk*(1)] + (u)[di*(0)+dj*(1)+dk*(-1)]) + 64*((u)[di*(0)+dj*(-1)+dk*(-1)] + (u)[di*(0)+dj*(1)+dk*(1)]) + 8*((u)[di*(0)+dj*(-1)+dk*(2)] + (u)[di*(0)+dj*(1)+dk*(-2)] + (u)[di*(0)+dj*(-2)+dk*(1)] + (u)[di*(0)+dj*(2)+dk*(-1)]) - 8*((u)[di*(0)+dj*(-1)+dk*(-2)] + (u)[di*(0)+dj*(1)+dk*(2)] + (u)[di*(0)+dj*(-2)+dk*(-1)] + (u)[di*(0)+dj*(2)+dk*(1)]) + (u)[di*(0)+dj*(-2)+dk*(-2)] - (u)[di*(0)+dj*(-2)+dk*(2)] - (u)[di*(0)+dj*(2)+dk*(-2)] + (u)[di*(0)+dj*(2)+dk*(2)]))
-#define DiffPlus1(u) (p1o1*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(1)+dj*(0)+dk*(0)]))
-#define DiffPlus2(u) (p1o1*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(0)+dj*(1)+dk*(0)]))
-#define DiffPlus3(u) (p1o1*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(0)+dj*(0)+dk*(1)]))
-#define DiffMinus1(u) (p1o1*((u)[di*(0)+dj*(0)+dk*(0)] - (u)[di*(-1)+dj*(0)+dk*(0)]))
-#define DiffMinus2(u) (p1o1*((u)[di*(0)+dj*(0)+dk*(0)] - (u)[di*(0)+dj*(-1)+dk*(0)]))
-#define DiffMinus3(u) (p1o1*((u)[di*(0)+dj*(0)+dk*(0)] - (u)[di*(0)+dj*(0)+dk*(-1)]))
-#define ShiftMinus1(u) (p1o1*(u)[di*(-1)+dj*(0)+dk*(0)])
-#define ShiftMinus2(u) (p1o1*(u)[di*(0)+dj*(-1)+dk*(0)])
-#define ShiftMinus3(u) (p1o1*(u)[di*(0)+dj*(0)+dk*(-1)])
-#define PDplus1(u) (p1odx*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(1)+dj*(0)+dk*(0)]))
-#define PDplus2(u) (p1ody*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(0)+dj*(1)+dk*(0)]))
-#define PDplus3(u) (p1odz*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(0)+dj*(0)+dk*(1)]))
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd1(u) ((-KRANC_GFOFFSET3D(u,-1,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1o2dx)
+#else
+#  define PDstandard2nd1(u) (PDstandard2nd1_impl(u,p1o2dx,cdj,cdk))
+static CCTK_REAL PDstandard2nd1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o2dx, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o2dx, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,-1,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1o2dx;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd2(u) ((-KRANC_GFOFFSET3D(u,0,-1,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1o2dy)
+#else
+#  define PDstandard2nd2(u) (PDstandard2nd2_impl(u,p1o2dy,cdj,cdk))
+static CCTK_REAL PDstandard2nd2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o2dy, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o2dy, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,-1,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1o2dy;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd3(u) ((-KRANC_GFOFFSET3D(u,0,0,-1) + KRANC_GFOFFSET3D(u,0,0,1))*p1o2dz)
+#else
+#  define PDstandard2nd3(u) (PDstandard2nd3_impl(u,p1o2dz,cdj,cdk))
+static CCTK_REAL PDstandard2nd3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o2dz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o2dz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,-1) + KRANC_GFOFFSET3D(u,0,0,1))*p1o2dz;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd11(u) ((-2*KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,-1,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1odx2)
+#else
+#  define PDstandard2nd11(u) (PDstandard2nd11_impl(u,p1odx2,cdj,cdk))
+static CCTK_REAL PDstandard2nd11_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odx2, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd11_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odx2, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-2*KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,-1,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1odx2;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd22(u) ((-2*KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,-1,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1ody2)
+#else
+#  define PDstandard2nd22(u) (PDstandard2nd22_impl(u,p1ody2,cdj,cdk))
+static CCTK_REAL PDstandard2nd22_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1ody2, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd22_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1ody2, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-2*KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,-1,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1ody2;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd33(u) ((-2*KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,0,-1) + KRANC_GFOFFSET3D(u,0,0,1))*p1odz2)
+#else
+#  define PDstandard2nd33(u) (PDstandard2nd33_impl(u,p1odz2,cdj,cdk))
+static CCTK_REAL PDstandard2nd33_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odz2, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd33_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odz2, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-2*KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,0,-1) + KRANC_GFOFFSET3D(u,0,0,1))*p1odz2;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd12(u) ((KRANC_GFOFFSET3D(u,-1,-1,0) - KRANC_GFOFFSET3D(u,-1,1,0) - KRANC_GFOFFSET3D(u,1,-1,0) + KRANC_GFOFFSET3D(u,1,1,0))*p1o4dxdy)
+#else
+#  define PDstandard2nd12(u) (PDstandard2nd12_impl(u,p1o4dxdy,cdj,cdk))
+static CCTK_REAL PDstandard2nd12_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dxdy, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd12_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dxdy, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,-1,-1,0) - KRANC_GFOFFSET3D(u,-1,1,0) - KRANC_GFOFFSET3D(u,1,-1,0) + KRANC_GFOFFSET3D(u,1,1,0))*p1o4dxdy;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd13(u) ((KRANC_GFOFFSET3D(u,-1,0,-1) - KRANC_GFOFFSET3D(u,-1,0,1) - KRANC_GFOFFSET3D(u,1,0,-1) + KRANC_GFOFFSET3D(u,1,0,1))*p1o4dxdz)
+#else
+#  define PDstandard2nd13(u) (PDstandard2nd13_impl(u,p1o4dxdz,cdj,cdk))
+static CCTK_REAL PDstandard2nd13_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dxdz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd13_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dxdz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,-1,0,-1) - KRANC_GFOFFSET3D(u,-1,0,1) - KRANC_GFOFFSET3D(u,1,0,-1) + KRANC_GFOFFSET3D(u,1,0,1))*p1o4dxdz;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd21(u) ((KRANC_GFOFFSET3D(u,-1,-1,0) - KRANC_GFOFFSET3D(u,-1,1,0) - KRANC_GFOFFSET3D(u,1,-1,0) + KRANC_GFOFFSET3D(u,1,1,0))*p1o4dxdy)
+#else
+#  define PDstandard2nd21(u) (PDstandard2nd21_impl(u,p1o4dxdy,cdj,cdk))
+static CCTK_REAL PDstandard2nd21_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dxdy, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd21_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dxdy, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,-1,-1,0) - KRANC_GFOFFSET3D(u,-1,1,0) - KRANC_GFOFFSET3D(u,1,-1,0) + KRANC_GFOFFSET3D(u,1,1,0))*p1o4dxdy;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd23(u) ((KRANC_GFOFFSET3D(u,0,-1,-1) - KRANC_GFOFFSET3D(u,0,-1,1) - KRANC_GFOFFSET3D(u,0,1,-1) + KRANC_GFOFFSET3D(u,0,1,1))*p1o4dydz)
+#else
+#  define PDstandard2nd23(u) (PDstandard2nd23_impl(u,p1o4dydz,cdj,cdk))
+static CCTK_REAL PDstandard2nd23_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dydz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd23_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dydz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,0,-1,-1) - KRANC_GFOFFSET3D(u,0,-1,1) - KRANC_GFOFFSET3D(u,0,1,-1) + KRANC_GFOFFSET3D(u,0,1,1))*p1o4dydz;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd31(u) ((KRANC_GFOFFSET3D(u,-1,0,-1) - KRANC_GFOFFSET3D(u,-1,0,1) - KRANC_GFOFFSET3D(u,1,0,-1) + KRANC_GFOFFSET3D(u,1,0,1))*p1o4dxdz)
+#else
+#  define PDstandard2nd31(u) (PDstandard2nd31_impl(u,p1o4dxdz,cdj,cdk))
+static CCTK_REAL PDstandard2nd31_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dxdz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd31_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dxdz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,-1,0,-1) - KRANC_GFOFFSET3D(u,-1,0,1) - KRANC_GFOFFSET3D(u,1,0,-1) + KRANC_GFOFFSET3D(u,1,0,1))*p1o4dxdz;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard2nd32(u) ((KRANC_GFOFFSET3D(u,0,-1,-1) - KRANC_GFOFFSET3D(u,0,-1,1) - KRANC_GFOFFSET3D(u,0,1,-1) + KRANC_GFOFFSET3D(u,0,1,1))*p1o4dydz)
+#else
+#  define PDstandard2nd32(u) (PDstandard2nd32_impl(u,p1o4dydz,cdj,cdk))
+static CCTK_REAL PDstandard2nd32_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dydz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard2nd32_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o4dydz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,0,-1,-1) - KRANC_GFOFFSET3D(u,0,-1,1) - KRANC_GFOFFSET3D(u,0,1,-1) + KRANC_GFOFFSET3D(u,0,1,1))*p1o4dydz;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th1(u) ((-8*KRANC_GFOFFSET3D(u,-1,0,0) + 8*KRANC_GFOFFSET3D(u,1,0,0) + KRANC_GFOFFSET3D(u,-2,0,0) - KRANC_GFOFFSET3D(u,2,0,0))*p1o12dx)
+#else
+#  define PDstandard4th1(u) (PDstandard4th1_impl(u,p1o12dx,cdj,cdk))
+static CCTK_REAL PDstandard4th1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o12dx, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o12dx, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-8*KRANC_GFOFFSET3D(u,-1,0,0) + 8*KRANC_GFOFFSET3D(u,1,0,0) + KRANC_GFOFFSET3D(u,-2,0,0) - KRANC_GFOFFSET3D(u,2,0,0))*p1o12dx;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th2(u) ((-8*KRANC_GFOFFSET3D(u,0,-1,0) + 8*KRANC_GFOFFSET3D(u,0,1,0) + KRANC_GFOFFSET3D(u,0,-2,0) - KRANC_GFOFFSET3D(u,0,2,0))*p1o12dy)
+#else
+#  define PDstandard4th2(u) (PDstandard4th2_impl(u,p1o12dy,cdj,cdk))
+static CCTK_REAL PDstandard4th2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o12dy, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o12dy, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-8*KRANC_GFOFFSET3D(u,0,-1,0) + 8*KRANC_GFOFFSET3D(u,0,1,0) + KRANC_GFOFFSET3D(u,0,-2,0) - KRANC_GFOFFSET3D(u,0,2,0))*p1o12dy;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th3(u) ((-8*KRANC_GFOFFSET3D(u,0,0,-1) + 8*KRANC_GFOFFSET3D(u,0,0,1) + KRANC_GFOFFSET3D(u,0,0,-2) - KRANC_GFOFFSET3D(u,0,0,2))*p1o12dz)
+#else
+#  define PDstandard4th3(u) (PDstandard4th3_impl(u,p1o12dz,cdj,cdk))
+static CCTK_REAL PDstandard4th3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o12dz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o12dz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-8*KRANC_GFOFFSET3D(u,0,0,-1) + 8*KRANC_GFOFFSET3D(u,0,0,1) + KRANC_GFOFFSET3D(u,0,0,-2) - KRANC_GFOFFSET3D(u,0,0,2))*p1o12dz;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th11(u) ((30*KRANC_GFOFFSET3D(u,0,0,0) - 16*(KRANC_GFOFFSET3D(u,-1,0,0) + KRANC_GFOFFSET3D(u,1,0,0)) + KRANC_GFOFFSET3D(u,-2,0,0) + KRANC_GFOFFSET3D(u,2,0,0))*pm1o12dx2)
+#else
+#  define PDstandard4th11(u) (PDstandard4th11_impl(u,pm1o12dx2,cdj,cdk))
+static CCTK_REAL PDstandard4th11_impl(CCTK_REAL const* restrict const u, CCTK_REAL const pm1o12dx2, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th11_impl(CCTK_REAL const* restrict const u, CCTK_REAL const pm1o12dx2, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (30*KRANC_GFOFFSET3D(u,0,0,0) - 16*(KRANC_GFOFFSET3D(u,-1,0,0) + KRANC_GFOFFSET3D(u,1,0,0)) + KRANC_GFOFFSET3D(u,-2,0,0) + KRANC_GFOFFSET3D(u,2,0,0))*pm1o12dx2;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th22(u) ((30*KRANC_GFOFFSET3D(u,0,0,0) - 16*(KRANC_GFOFFSET3D(u,0,-1,0) + KRANC_GFOFFSET3D(u,0,1,0)) + KRANC_GFOFFSET3D(u,0,-2,0) + KRANC_GFOFFSET3D(u,0,2,0))*pm1o12dy2)
+#else
+#  define PDstandard4th22(u) (PDstandard4th22_impl(u,pm1o12dy2,cdj,cdk))
+static CCTK_REAL PDstandard4th22_impl(CCTK_REAL const* restrict const u, CCTK_REAL const pm1o12dy2, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th22_impl(CCTK_REAL const* restrict const u, CCTK_REAL const pm1o12dy2, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (30*KRANC_GFOFFSET3D(u,0,0,0) - 16*(KRANC_GFOFFSET3D(u,0,-1,0) + KRANC_GFOFFSET3D(u,0,1,0)) + KRANC_GFOFFSET3D(u,0,-2,0) + KRANC_GFOFFSET3D(u,0,2,0))*pm1o12dy2;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th33(u) ((30*KRANC_GFOFFSET3D(u,0,0,0) - 16*(KRANC_GFOFFSET3D(u,0,0,-1) + KRANC_GFOFFSET3D(u,0,0,1)) + KRANC_GFOFFSET3D(u,0,0,-2) + KRANC_GFOFFSET3D(u,0,0,2))*pm1o12dz2)
+#else
+#  define PDstandard4th33(u) (PDstandard4th33_impl(u,pm1o12dz2,cdj,cdk))
+static CCTK_REAL PDstandard4th33_impl(CCTK_REAL const* restrict const u, CCTK_REAL const pm1o12dz2, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th33_impl(CCTK_REAL const* restrict const u, CCTK_REAL const pm1o12dz2, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (30*KRANC_GFOFFSET3D(u,0,0,0) - 16*(KRANC_GFOFFSET3D(u,0,0,-1) + KRANC_GFOFFSET3D(u,0,0,1)) + KRANC_GFOFFSET3D(u,0,0,-2) + KRANC_GFOFFSET3D(u,0,0,2))*pm1o12dz2;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th12(u) ((-64*(KRANC_GFOFFSET3D(u,-1,1,0) + KRANC_GFOFFSET3D(u,1,-1,0)) + 64*(KRANC_GFOFFSET3D(u,-1,-1,0) + KRANC_GFOFFSET3D(u,1,1,0)) + 8*(KRANC_GFOFFSET3D(u,-1,2,0) + KRANC_GFOFFSET3D(u,1,-2,0) + KRANC_GFOFFSET3D(u,-2,1,0) + KRANC_GFOFFSET3D(u,2,-1,0)) - 8*(KRANC_GFOFFSET3D(u,-1,-2,0) + KRANC_GFOFFSET3D(u,1,2,0) + KRANC_GFOFFSET3D(u,-2,-1,0) + KRANC_GFOFFSET3D(u,2,1,0)) + KRANC_GFOFFSET3D(u,-2,-2,0) - KRANC_GFOFFSET3D(u,-2,2,0) - KRANC_GFOFFSET3D(u,2,-2,0) + KRANC_GFOFFSET3D(u,2,2,0))*p1o144dxdy)
+#else
+#  define PDstandard4th12(u) (PDstandard4th12_impl(u,p1o144dxdy,cdj,cdk))
+static CCTK_REAL PDstandard4th12_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dxdy, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th12_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dxdy, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-64*(KRANC_GFOFFSET3D(u,-1,1,0) + KRANC_GFOFFSET3D(u,1,-1,0)) + 64*(KRANC_GFOFFSET3D(u,-1,-1,0) + KRANC_GFOFFSET3D(u,1,1,0)) + 8*(KRANC_GFOFFSET3D(u,-1,2,0) + KRANC_GFOFFSET3D(u,1,-2,0) + KRANC_GFOFFSET3D(u,-2,1,0) + KRANC_GFOFFSET3D(u,2,-1,0)) - 8*(KRANC_GFOFFSET3D(u,-1,-2,0) + KRANC_GFOFFSET3D(u,1,2,0) + KRANC_GFOFFSET3D(u,-2,-1,0) + KRANC_GFOFFSET3D(u,2,1,0)) + KRANC_GFOFFSET3D(u,-2,-2,0) - KRANC_GFOFFSET3D(u,-2,2,0) - KRANC_GFOFFSET3D(u,2,-2,0) + KRANC_GFOFFSET3D(u,2,2,0))*p1o144dxdy;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th13(u) ((-64*(KRANC_GFOFFSET3D(u,-1,0,1) + KRANC_GFOFFSET3D(u,1,0,-1)) + 64*(KRANC_GFOFFSET3D(u,-1,0,-1) + KRANC_GFOFFSET3D(u,1,0,1)) + 8*(KRANC_GFOFFSET3D(u,-1,0,2) + KRANC_GFOFFSET3D(u,1,0,-2) + KRANC_GFOFFSET3D(u,-2,0,1) + KRANC_GFOFFSET3D(u,2,0,-1)) - 8*(KRANC_GFOFFSET3D(u,-1,0,-2) + KRANC_GFOFFSET3D(u,1,0,2) + KRANC_GFOFFSET3D(u,-2,0,-1) + KRANC_GFOFFSET3D(u,2,0,1)) + KRANC_GFOFFSET3D(u,-2,0,-2) - KRANC_GFOFFSET3D(u,-2,0,2) - KRANC_GFOFFSET3D(u,2,0,-2) + KRANC_GFOFFSET3D(u,2,0,2))*p1o144dxdz)
+#else
+#  define PDstandard4th13(u) (PDstandard4th13_impl(u,p1o144dxdz,cdj,cdk))
+static CCTK_REAL PDstandard4th13_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dxdz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th13_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dxdz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-64*(KRANC_GFOFFSET3D(u,-1,0,1) + KRANC_GFOFFSET3D(u,1,0,-1)) + 64*(KRANC_GFOFFSET3D(u,-1,0,-1) + KRANC_GFOFFSET3D(u,1,0,1)) + 8*(KRANC_GFOFFSET3D(u,-1,0,2) + KRANC_GFOFFSET3D(u,1,0,-2) + KRANC_GFOFFSET3D(u,-2,0,1) + KRANC_GFOFFSET3D(u,2,0,-1)) - 8*(KRANC_GFOFFSET3D(u,-1,0,-2) + KRANC_GFOFFSET3D(u,1,0,2) + KRANC_GFOFFSET3D(u,-2,0,-1) + KRANC_GFOFFSET3D(u,2,0,1)) + KRANC_GFOFFSET3D(u,-2,0,-2) - KRANC_GFOFFSET3D(u,-2,0,2) - KRANC_GFOFFSET3D(u,2,0,-2) + KRANC_GFOFFSET3D(u,2,0,2))*p1o144dxdz;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th21(u) ((-64*(KRANC_GFOFFSET3D(u,-1,1,0) + KRANC_GFOFFSET3D(u,1,-1,0)) + 64*(KRANC_GFOFFSET3D(u,-1,-1,0) + KRANC_GFOFFSET3D(u,1,1,0)) + 8*(KRANC_GFOFFSET3D(u,-1,2,0) + KRANC_GFOFFSET3D(u,1,-2,0) + KRANC_GFOFFSET3D(u,-2,1,0) + KRANC_GFOFFSET3D(u,2,-1,0)) - 8*(KRANC_GFOFFSET3D(u,-1,-2,0) + KRANC_GFOFFSET3D(u,1,2,0) + KRANC_GFOFFSET3D(u,-2,-1,0) + KRANC_GFOFFSET3D(u,2,1,0)) + KRANC_GFOFFSET3D(u,-2,-2,0) - KRANC_GFOFFSET3D(u,-2,2,0) - KRANC_GFOFFSET3D(u,2,-2,0) + KRANC_GFOFFSET3D(u,2,2,0))*p1o144dxdy)
+#else
+#  define PDstandard4th21(u) (PDstandard4th21_impl(u,p1o144dxdy,cdj,cdk))
+static CCTK_REAL PDstandard4th21_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dxdy, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th21_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dxdy, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-64*(KRANC_GFOFFSET3D(u,-1,1,0) + KRANC_GFOFFSET3D(u,1,-1,0)) + 64*(KRANC_GFOFFSET3D(u,-1,-1,0) + KRANC_GFOFFSET3D(u,1,1,0)) + 8*(KRANC_GFOFFSET3D(u,-1,2,0) + KRANC_GFOFFSET3D(u,1,-2,0) + KRANC_GFOFFSET3D(u,-2,1,0) + KRANC_GFOFFSET3D(u,2,-1,0)) - 8*(KRANC_GFOFFSET3D(u,-1,-2,0) + KRANC_GFOFFSET3D(u,1,2,0) + KRANC_GFOFFSET3D(u,-2,-1,0) + KRANC_GFOFFSET3D(u,2,1,0)) + KRANC_GFOFFSET3D(u,-2,-2,0) - KRANC_GFOFFSET3D(u,-2,2,0) - KRANC_GFOFFSET3D(u,2,-2,0) + KRANC_GFOFFSET3D(u,2,2,0))*p1o144dxdy;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th23(u) ((-64*(KRANC_GFOFFSET3D(u,0,-1,1) + KRANC_GFOFFSET3D(u,0,1,-1)) + 64*(KRANC_GFOFFSET3D(u,0,-1,-1) + KRANC_GFOFFSET3D(u,0,1,1)) + 8*(KRANC_GFOFFSET3D(u,0,-1,2) + KRANC_GFOFFSET3D(u,0,1,-2) + KRANC_GFOFFSET3D(u,0,-2,1) + KRANC_GFOFFSET3D(u,0,2,-1)) - 8*(KRANC_GFOFFSET3D(u,0,-1,-2) + KRANC_GFOFFSET3D(u,0,1,2) + KRANC_GFOFFSET3D(u,0,-2,-1) + KRANC_GFOFFSET3D(u,0,2,1)) + KRANC_GFOFFSET3D(u,0,-2,-2) - KRANC_GFOFFSET3D(u,0,-2,2) - KRANC_GFOFFSET3D(u,0,2,-2) + KRANC_GFOFFSET3D(u,0,2,2))*p1o144dydz)
+#else
+#  define PDstandard4th23(u) (PDstandard4th23_impl(u,p1o144dydz,cdj,cdk))
+static CCTK_REAL PDstandard4th23_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dydz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th23_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dydz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-64*(KRANC_GFOFFSET3D(u,0,-1,1) + KRANC_GFOFFSET3D(u,0,1,-1)) + 64*(KRANC_GFOFFSET3D(u,0,-1,-1) + KRANC_GFOFFSET3D(u,0,1,1)) + 8*(KRANC_GFOFFSET3D(u,0,-1,2) + KRANC_GFOFFSET3D(u,0,1,-2) + KRANC_GFOFFSET3D(u,0,-2,1) + KRANC_GFOFFSET3D(u,0,2,-1)) - 8*(KRANC_GFOFFSET3D(u,0,-1,-2) + KRANC_GFOFFSET3D(u,0,1,2) + KRANC_GFOFFSET3D(u,0,-2,-1) + KRANC_GFOFFSET3D(u,0,2,1)) + KRANC_GFOFFSET3D(u,0,-2,-2) - KRANC_GFOFFSET3D(u,0,-2,2) - KRANC_GFOFFSET3D(u,0,2,-2) + KRANC_GFOFFSET3D(u,0,2,2))*p1o144dydz;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th31(u) ((-64*(KRANC_GFOFFSET3D(u,-1,0,1) + KRANC_GFOFFSET3D(u,1,0,-1)) + 64*(KRANC_GFOFFSET3D(u,-1,0,-1) + KRANC_GFOFFSET3D(u,1,0,1)) + 8*(KRANC_GFOFFSET3D(u,-1,0,2) + KRANC_GFOFFSET3D(u,1,0,-2) + KRANC_GFOFFSET3D(u,-2,0,1) + KRANC_GFOFFSET3D(u,2,0,-1)) - 8*(KRANC_GFOFFSET3D(u,-1,0,-2) + KRANC_GFOFFSET3D(u,1,0,2) + KRANC_GFOFFSET3D(u,-2,0,-1) + KRANC_GFOFFSET3D(u,2,0,1)) + KRANC_GFOFFSET3D(u,-2,0,-2) - KRANC_GFOFFSET3D(u,-2,0,2) - KRANC_GFOFFSET3D(u,2,0,-2) + KRANC_GFOFFSET3D(u,2,0,2))*p1o144dxdz)
+#else
+#  define PDstandard4th31(u) (PDstandard4th31_impl(u,p1o144dxdz,cdj,cdk))
+static CCTK_REAL PDstandard4th31_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dxdz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th31_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dxdz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-64*(KRANC_GFOFFSET3D(u,-1,0,1) + KRANC_GFOFFSET3D(u,1,0,-1)) + 64*(KRANC_GFOFFSET3D(u,-1,0,-1) + KRANC_GFOFFSET3D(u,1,0,1)) + 8*(KRANC_GFOFFSET3D(u,-1,0,2) + KRANC_GFOFFSET3D(u,1,0,-2) + KRANC_GFOFFSET3D(u,-2,0,1) + KRANC_GFOFFSET3D(u,2,0,-1)) - 8*(KRANC_GFOFFSET3D(u,-1,0,-2) + KRANC_GFOFFSET3D(u,1,0,2) + KRANC_GFOFFSET3D(u,-2,0,-1) + KRANC_GFOFFSET3D(u,2,0,1)) + KRANC_GFOFFSET3D(u,-2,0,-2) - KRANC_GFOFFSET3D(u,-2,0,2) - KRANC_GFOFFSET3D(u,2,0,-2) + KRANC_GFOFFSET3D(u,2,0,2))*p1o144dxdz;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDstandard4th32(u) ((-64*(KRANC_GFOFFSET3D(u,0,-1,1) + KRANC_GFOFFSET3D(u,0,1,-1)) + 64*(KRANC_GFOFFSET3D(u,0,-1,-1) + KRANC_GFOFFSET3D(u,0,1,1)) + 8*(KRANC_GFOFFSET3D(u,0,-1,2) + KRANC_GFOFFSET3D(u,0,1,-2) + KRANC_GFOFFSET3D(u,0,-2,1) + KRANC_GFOFFSET3D(u,0,2,-1)) - 8*(KRANC_GFOFFSET3D(u,0,-1,-2) + KRANC_GFOFFSET3D(u,0,1,2) + KRANC_GFOFFSET3D(u,0,-2,-1) + KRANC_GFOFFSET3D(u,0,2,1)) + KRANC_GFOFFSET3D(u,0,-2,-2) - KRANC_GFOFFSET3D(u,0,-2,2) - KRANC_GFOFFSET3D(u,0,2,-2) + KRANC_GFOFFSET3D(u,0,2,2))*p1o144dydz)
+#else
+#  define PDstandard4th32(u) (PDstandard4th32_impl(u,p1o144dydz,cdj,cdk))
+static CCTK_REAL PDstandard4th32_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dydz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDstandard4th32_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o144dydz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-64*(KRANC_GFOFFSET3D(u,0,-1,1) + KRANC_GFOFFSET3D(u,0,1,-1)) + 64*(KRANC_GFOFFSET3D(u,0,-1,-1) + KRANC_GFOFFSET3D(u,0,1,1)) + 8*(KRANC_GFOFFSET3D(u,0,-1,2) + KRANC_GFOFFSET3D(u,0,1,-2) + KRANC_GFOFFSET3D(u,0,-2,1) + KRANC_GFOFFSET3D(u,0,2,-1)) - 8*(KRANC_GFOFFSET3D(u,0,-1,-2) + KRANC_GFOFFSET3D(u,0,1,2) + KRANC_GFOFFSET3D(u,0,-2,-1) + KRANC_GFOFFSET3D(u,0,2,1)) + KRANC_GFOFFSET3D(u,0,-2,-2) - KRANC_GFOFFSET3D(u,0,-2,2) - KRANC_GFOFFSET3D(u,0,2,-2) + KRANC_GFOFFSET3D(u,0,2,2))*p1o144dydz;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffPlus1(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1o1)
+#else
+#  define DiffPlus1(u) (DiffPlus1_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffPlus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffPlus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffPlus2(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1o1)
+#else
+#  define DiffPlus2(u) (DiffPlus2_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffPlus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffPlus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffPlus3(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,0,1))*p1o1)
+#else
+#  define DiffPlus3(u) (DiffPlus3_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffPlus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffPlus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,0,1))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffMinus1(u) ((KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,-1,0,0))*p1o1)
+#else
+#  define DiffMinus1(u) (DiffMinus1_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffMinus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffMinus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,-1,0,0))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffMinus2(u) ((KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,0,-1,0))*p1o1)
+#else
+#  define DiffMinus2(u) (DiffMinus2_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffMinus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffMinus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,0,-1,0))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffMinus3(u) ((KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,0,0,-1))*p1o1)
+#else
+#  define DiffMinus3(u) (DiffMinus3_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffMinus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffMinus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,0,0,-1))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define ShiftMinus1(u) (KRANC_GFOFFSET3D(u,-1,0,0)*p1o1)
+#else
+#  define ShiftMinus1(u) (ShiftMinus1_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL ShiftMinus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL ShiftMinus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return KRANC_GFOFFSET3D(u,-1,0,0)*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define ShiftMinus2(u) (KRANC_GFOFFSET3D(u,0,-1,0)*p1o1)
+#else
+#  define ShiftMinus2(u) (ShiftMinus2_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL ShiftMinus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL ShiftMinus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return KRANC_GFOFFSET3D(u,0,-1,0)*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define ShiftMinus3(u) (KRANC_GFOFFSET3D(u,0,0,-1)*p1o1)
+#else
+#  define ShiftMinus3(u) (ShiftMinus3_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL ShiftMinus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL ShiftMinus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return KRANC_GFOFFSET3D(u,0,0,-1)*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDplus1(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1odx)
+#else
+#  define PDplus1(u) (PDplus1_impl(u,p1odx,cdj,cdk))
+static CCTK_REAL PDplus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odx, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDplus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odx, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1odx;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDplus2(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1ody)
+#else
+#  define PDplus2(u) (PDplus2_impl(u,p1ody,cdj,cdk))
+static CCTK_REAL PDplus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1ody, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDplus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1ody, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1ody;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDplus3(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,0,1))*p1odz)
+#else
+#  define PDplus3(u) (PDplus3_impl(u,p1odz,cdj,cdk))
+static CCTK_REAL PDplus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDplus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,0,1))*p1odz;
+}
+#endif
+
diff --git a/Examples/EM/src/EM_constraints.cc b/Examples/EM/src/EM_constraints.cc
index d26cb66..dc40783 100644
--- a/Examples/EM/src/EM_constraints.cc
+++ b/Examples/EM/src/EM_constraints.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -20,43 +22,40 @@
 #define SQR(x) ((x) * (x))
 #define CUB(x) ((x) * (x) * (x))
 
-static void EM_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[])
+extern "C" void EM_constraints_SelectBCs(CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  CCTK_INT ierr = 0;
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "EM::constraints","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for EM::constraints.");
+  return;
+}
+
+static void EM_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 imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
   
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering EM_constraints_Body");
-  }
-  
-  if (cctk_iteration % EM_constraints_calc_every != EM_constraints_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EM::B_group","EM::constraints","EM::El_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "EM_constraints", 3, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "EM_constraints", 1, 1, 1);
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -93,48 +92,54 @@ static void EM_constraints_Body(cGH const * restrict const cctkGH, int const dir
   CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy));
   CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz));
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (EM_constraints,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL B1L = B1[index];
-        CCTK_REAL B2L = B2[index];
-        CCTK_REAL B3L = B3[index];
-        CCTK_REAL El1L = El1[index];
-        CCTK_REAL El2L = El2[index];
-        CCTK_REAL El3L = El3[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const PDstandard2nd1B1 = PDstandard2nd1(&B1[index]);
-        CCTK_REAL const PDstandard2nd2B2 = PDstandard2nd2(&B2[index]);
-        CCTK_REAL const PDstandard2nd3B3 = PDstandard2nd3(&B3[index]);
-        CCTK_REAL const PDstandard2nd1El1 = PDstandard2nd1(&El1[index]);
-        CCTK_REAL const PDstandard2nd2El2 = PDstandard2nd2(&El2[index]);
-        CCTK_REAL const PDstandard2nd3El3 = PDstandard2nd3(&El3[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL CElL = PDstandard2nd1El1 + PDstandard2nd2El2 + 
-          PDstandard2nd3El3;
-        
-        CCTK_REAL CBL = PDstandard2nd1B1 + PDstandard2nd2B2 + 
-          PDstandard2nd3B3;
-        
-        /* Copy local copies back to grid functions */
-        CB[index] = CBL;
-        CEl[index] = CElL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL B1L = B1[index];
+    CCTK_REAL B2L = B2[index];
+    CCTK_REAL B3L = B3[index];
+    CCTK_REAL El1L = El1[index];
+    CCTK_REAL El2L = El2[index];
+    CCTK_REAL El3L = El3[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDstandard2nd1B1 = PDstandard2nd1(&B1[index]);
+    CCTK_REAL const PDstandard2nd2B2 = PDstandard2nd2(&B2[index]);
+    CCTK_REAL const PDstandard2nd3B3 = PDstandard2nd3(&B3[index]);
+    CCTK_REAL const PDstandard2nd1El1 = PDstandard2nd1(&El1[index]);
+    CCTK_REAL const PDstandard2nd2El2 = PDstandard2nd2(&El2[index]);
+    CCTK_REAL const PDstandard2nd3El3 = PDstandard2nd3(&El3[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL CElL = PDstandard2nd1El1 + PDstandard2nd2El2 + 
+      PDstandard2nd3El3;
+    
+    CCTK_REAL CBL = PDstandard2nd1B1 + PDstandard2nd2B2 + 
+      PDstandard2nd3B3;
+    
+    /* Copy local copies back to grid functions */
+    CB[index] = CBL;
+    CEl[index] = CElL;
   }
+  CCTK_ENDLOOP3 (EM_constraints);
 }
 
 extern "C" void EM_constraints(CCTK_ARGUMENTS)
@@ -142,5 +147,26 @@ extern "C" void EM_constraints(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
-  GenericFD_LoopOverEverything(cctkGH, &EM_constraints_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering EM_constraints_Body");
+  }
+  
+  if (cctk_iteration % EM_constraints_calc_every != EM_constraints_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EM::B_group","EM::constraints","EM::El_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "EM_constraints", 3, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "EM_constraints", 1, 1, 1);
+  
+  GenericFD_LoopOverInterior(cctkGH, &EM_constraints_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving EM_constraints_Body");
+  }
 }
diff --git a/Examples/EM/src/EM_energy.cc b/Examples/EM/src/EM_energy.cc
index 2318bd8..b8e1a17 100644
--- a/Examples/EM/src/EM_energy.cc
+++ b/Examples/EM/src/EM_energy.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -20,42 +22,28 @@
 #define SQR(x) ((x) * (x))
 #define CUB(x) ((x) * (x) * (x))
 
-static void EM_energy_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 EM_energy_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 imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering EM_energy_Body");
-  }
-  
-  if (cctk_iteration % EM_energy_calc_every != EM_energy_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EM::B_group","EM::El_group","EM::endens"};
-  GenericFD_AssertGroupStorage(cctkGH, "EM_energy", 3, groups);
-  
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -92,38 +80,44 @@ static void EM_energy_Body(cGH const * restrict const cctkGH, int const dir, int
   CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy));
   CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz));
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (EM_energy,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL B1L = B1[index];
-        CCTK_REAL B2L = B2[index];
-        CCTK_REAL B3L = B3[index];
-        CCTK_REAL El1L = El1[index];
-        CCTK_REAL El2L = El2[index];
-        CCTK_REAL El3L = El3[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL rhoL = 0.5*(SQR(B1L) + SQR(B2L) + SQR(B3L) + SQR(El1L) + 
-          SQR(El2L) + SQR(El3L));
-        
-        /* Copy local copies back to grid functions */
-        rho[index] = rhoL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL B1L = B1[index];
+    CCTK_REAL B2L = B2[index];
+    CCTK_REAL B3L = B3[index];
+    CCTK_REAL El1L = El1[index];
+    CCTK_REAL El2L = El2[index];
+    CCTK_REAL El3L = El3[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL rhoL = 0.5*(SQR(B1L) + SQR(B2L) + SQR(B3L) + 
+      SQR(El1L) + SQR(El2L) + SQR(El3L));
+    
+    /* Copy local copies back to grid functions */
+    rho[index] = rhoL;
   }
+  CCTK_ENDLOOP3 (EM_energy);
 }
 
 extern "C" void EM_energy(CCTK_ARGUMENTS)
@@ -131,5 +125,25 @@ extern "C" void EM_energy(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering EM_energy_Body");
+  }
+  
+  if (cctk_iteration % EM_energy_calc_every != EM_energy_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EM::B_group","EM::El_group","EM::endens"};
+  GenericFD_AssertGroupStorage(cctkGH, "EM_energy", 3, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &EM_energy_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving EM_energy_Body");
+  }
 }
diff --git a/Examples/EM/src/EM_evol.cc b/Examples/EM/src/EM_evol.cc
index 311d5ca..ad67af3 100644
--- a/Examples/EM/src/EM_evol.cc
+++ b/Examples/EM/src/EM_evol.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -20,43 +22,43 @@
 #define SQR(x) ((x) * (x))
 #define CUB(x) ((x) * (x) * (x))
 
-static void EM_evol_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[])
+extern "C" void EM_evol_SelectBCs(CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  CCTK_INT ierr = 0;
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "EM::B_grouprhs","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for EM::B_grouprhs.");
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "EM::El_grouprhs","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for EM::El_grouprhs.");
+  return;
+}
+
+static void EM_evol_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 imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
   
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering EM_evol_Body");
-  }
-  
-  if (cctk_iteration % EM_evol_calc_every != EM_evol_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EM::B_group","EM::B_grouprhs","EM::El_group","EM::El_grouprhs"};
-  GenericFD_AssertGroupStorage(cctkGH, "EM_evol", 4, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "EM_evol", 1, 1, 1);
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -93,64 +95,70 @@ static void EM_evol_Body(cGH const * restrict const cctkGH, int const dir, int c
   CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy));
   CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz));
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (EM_evol,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL B1L = B1[index];
-        CCTK_REAL B2L = B2[index];
-        CCTK_REAL B3L = B3[index];
-        CCTK_REAL El1L = El1[index];
-        CCTK_REAL El2L = El2[index];
-        CCTK_REAL El3L = El3[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const PDstandard2nd2B1 = PDstandard2nd2(&B1[index]);
-        CCTK_REAL const PDstandard2nd3B1 = PDstandard2nd3(&B1[index]);
-        CCTK_REAL const PDstandard2nd1B2 = PDstandard2nd1(&B2[index]);
-        CCTK_REAL const PDstandard2nd3B2 = PDstandard2nd3(&B2[index]);
-        CCTK_REAL const PDstandard2nd1B3 = PDstandard2nd1(&B3[index]);
-        CCTK_REAL const PDstandard2nd2B3 = PDstandard2nd2(&B3[index]);
-        CCTK_REAL const PDstandard2nd2El1 = PDstandard2nd2(&El1[index]);
-        CCTK_REAL const PDstandard2nd3El1 = PDstandard2nd3(&El1[index]);
-        CCTK_REAL const PDstandard2nd1El2 = PDstandard2nd1(&El2[index]);
-        CCTK_REAL const PDstandard2nd3El2 = PDstandard2nd3(&El2[index]);
-        CCTK_REAL const PDstandard2nd1El3 = PDstandard2nd1(&El3[index]);
-        CCTK_REAL const PDstandard2nd2El3 = PDstandard2nd2(&El3[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL El1rhsL = PDstandard2nd2B3 - PDstandard2nd3B2;
-        
-        CCTK_REAL El2rhsL = -PDstandard2nd1B3 + PDstandard2nd3B1;
-        
-        CCTK_REAL El3rhsL = PDstandard2nd1B2 - PDstandard2nd2B1;
-        
-        CCTK_REAL B1rhsL = -PDstandard2nd2El3 + PDstandard2nd3El2;
-        
-        CCTK_REAL B2rhsL = PDstandard2nd1El3 - PDstandard2nd3El1;
-        
-        CCTK_REAL B3rhsL = -PDstandard2nd1El2 + PDstandard2nd2El1;
-        
-        /* Copy local copies back to grid functions */
-        B1rhs[index] = B1rhsL;
-        B2rhs[index] = B2rhsL;
-        B3rhs[index] = B3rhsL;
-        El1rhs[index] = El1rhsL;
-        El2rhs[index] = El2rhsL;
-        El3rhs[index] = El3rhsL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL B1L = B1[index];
+    CCTK_REAL B2L = B2[index];
+    CCTK_REAL B3L = B3[index];
+    CCTK_REAL El1L = El1[index];
+    CCTK_REAL El2L = El2[index];
+    CCTK_REAL El3L = El3[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDstandard2nd2B1 = PDstandard2nd2(&B1[index]);
+    CCTK_REAL const PDstandard2nd3B1 = PDstandard2nd3(&B1[index]);
+    CCTK_REAL const PDstandard2nd1B2 = PDstandard2nd1(&B2[index]);
+    CCTK_REAL const PDstandard2nd3B2 = PDstandard2nd3(&B2[index]);
+    CCTK_REAL const PDstandard2nd1B3 = PDstandard2nd1(&B3[index]);
+    CCTK_REAL const PDstandard2nd2B3 = PDstandard2nd2(&B3[index]);
+    CCTK_REAL const PDstandard2nd2El1 = PDstandard2nd2(&El1[index]);
+    CCTK_REAL const PDstandard2nd3El1 = PDstandard2nd3(&El1[index]);
+    CCTK_REAL const PDstandard2nd1El2 = PDstandard2nd1(&El2[index]);
+    CCTK_REAL const PDstandard2nd3El2 = PDstandard2nd3(&El2[index]);
+    CCTK_REAL const PDstandard2nd1El3 = PDstandard2nd1(&El3[index]);
+    CCTK_REAL const PDstandard2nd2El3 = PDstandard2nd2(&El3[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL El1rhsL = PDstandard2nd2B3 - PDstandard2nd3B2;
+    
+    CCTK_REAL El2rhsL = -PDstandard2nd1B3 + PDstandard2nd3B1;
+    
+    CCTK_REAL El3rhsL = PDstandard2nd1B2 - PDstandard2nd2B1;
+    
+    CCTK_REAL B1rhsL = -PDstandard2nd2El3 + PDstandard2nd3El2;
+    
+    CCTK_REAL B2rhsL = PDstandard2nd1El3 - PDstandard2nd3El1;
+    
+    CCTK_REAL B3rhsL = -PDstandard2nd1El2 + PDstandard2nd2El1;
+    
+    /* Copy local copies back to grid functions */
+    B1rhs[index] = B1rhsL;
+    B2rhs[index] = B2rhsL;
+    B3rhs[index] = B3rhsL;
+    El1rhs[index] = El1rhsL;
+    El2rhs[index] = El2rhsL;
+    El3rhs[index] = El3rhsL;
   }
+  CCTK_ENDLOOP3 (EM_evol);
 }
 
 extern "C" void EM_evol(CCTK_ARGUMENTS)
@@ -158,5 +166,26 @@ extern "C" void EM_evol(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
-  GenericFD_LoopOverEverything(cctkGH, &EM_evol_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering EM_evol_Body");
+  }
+  
+  if (cctk_iteration % EM_evol_calc_every != EM_evol_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EM::B_group","EM::B_grouprhs","EM::El_group","EM::El_grouprhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "EM_evol", 4, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "EM_evol", 1, 1, 1);
+  
+  GenericFD_LoopOverInterior(cctkGH, &EM_evol_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving EM_evol_Body");
+  }
 }
diff --git a/Examples/EM/src/EM_initial.cc b/Examples/EM/src/EM_initial.cc
index 9b2f83e..8d0e0b0 100644
--- a/Examples/EM/src/EM_initial.cc
+++ b/Examples/EM/src/EM_initial.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -20,42 +22,28 @@
 #define SQR(x) ((x) * (x))
 #define CUB(x) ((x) * (x) * (x))
 
-static void EM_initial_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 EM_initial_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 imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering EM_initial_Body");
-  }
-  
-  if (cctk_iteration % EM_initial_calc_every != EM_initial_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EM::B_group","EM::El_group","grid::coordinates"};
-  GenericFD_AssertGroupStorage(cctkGH, "EM_initial", 3, groups);
-  
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -92,50 +80,56 @@ static void EM_initial_Body(cGH const * restrict const cctkGH, int const dir, in
   CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy));
   CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz));
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (EM_initial,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL xL = x[index];
-        CCTK_REAL yL = y[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL El1L = Cos(2*Pi*(xL + yL))*ToReal(sigma);
-        
-        CCTK_REAL El2L = Cos(2*Pi*xL)*(-1 + ToReal(sigma)) - Cos(2*Pi*(xL + 
-          yL))*ToReal(sigma);
-        
-        CCTK_REAL El3L = 0;
-        
-        CCTK_REAL B1L = 0;
-        
-        CCTK_REAL B2L = 0;
-        
-        CCTK_REAL B3L = -(Cos(2*Pi*xL)*(-1 + ToReal(sigma))) + Cos(2*Pi*(xL + 
-          yL))*ToReal(sigma);
-        
-        /* Copy local copies back to grid functions */
-        B1[index] = B1L;
-        B2[index] = B2L;
-        B3[index] = B3L;
-        El1[index] = El1L;
-        El2[index] = El2L;
-        El3[index] = El3L;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL xL = x[index];
+    CCTK_REAL yL = y[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL El1L = Cos(2*(xL + yL)*Pi)*ToReal(sigma);
+    
+    CCTK_REAL El2L = Cos(2*xL*Pi)*(-1 + ToReal(sigma)) - Cos(2*(xL + 
+      yL)*Pi)*ToReal(sigma);
+    
+    CCTK_REAL El3L = 0;
+    
+    CCTK_REAL B1L = 0;
+    
+    CCTK_REAL B2L = 0;
+    
+    CCTK_REAL B3L = -(Cos(2*xL*Pi)*(-1 + ToReal(sigma))) + Cos(2*(xL + 
+      yL)*Pi)*ToReal(sigma);
+    
+    /* Copy local copies back to grid functions */
+    B1[index] = B1L;
+    B2[index] = B2L;
+    B3[index] = B3L;
+    El1[index] = El1L;
+    El2[index] = El2L;
+    El3[index] = El3L;
   }
+  CCTK_ENDLOOP3 (EM_initial);
 }
 
 extern "C" void EM_initial(CCTK_ARGUMENTS)
@@ -143,5 +137,25 @@ extern "C" void EM_initial(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering EM_initial_Body");
+  }
+  
+  if (cctk_iteration % EM_initial_calc_every != EM_initial_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EM::B_group","EM::El_group","grid::coordinates"};
+  GenericFD_AssertGroupStorage(cctkGH, "EM_initial", 3, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &EM_initial_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving EM_initial_Body");
+  }
 }
diff --git a/Examples/EM/src/RegisterMoL.cc b/Examples/EM/src/RegisterMoL.cc
index 34583ae..02f967a 100644
--- a/Examples/EM/src/RegisterMoL.cc
+++ b/Examples/EM/src/RegisterMoL.cc
@@ -18,5 +18,7 @@ extern "C" void EM_RegisterVars(CCTK_ARGUMENTS)
   ierr += MoLRegisterEvolved(CCTK_VarIndex("EM::El1"),  CCTK_VarIndex("EM::El1rhs"));
   ierr += MoLRegisterEvolved(CCTK_VarIndex("EM::El2"),  CCTK_VarIndex("EM::El2rhs"));
   ierr += MoLRegisterEvolved(CCTK_VarIndex("EM::El3"),  CCTK_VarIndex("EM::El3rhs"));
+  
+  /* Register all the evolved Array functions with MoL */
   return;
 }