CarpetLib: Add preliminary support for DGFE operators

8 files changed, 574 insertions, 4 deletions

diff --git a/Carpet/CarpetLib/configuration.ccl b/Carpet/CarpetLib/configuration.ccl
index 313afc304..7e9f0f8bd 100644
--- a/Carpet/CarpetLib/configuration.ccl
+++ b/Carpet/CarpetLib/configuration.ccl
@@ -6,8 +6,13 @@ PROVIDES CarpetLib
   LANG
 }
 
+REQUIRES Vectors
+
 OPTIONAL LoopControl
 {
 }
 
-REQUIRES Vectors
+# DGFE
+OPTIONAL HRSCCore
+{
+}
diff --git a/Carpet/CarpetLib/interface.ccl b/Carpet/CarpetLib/interface.ccl
index e2184abeb..66ebc9f5d 100644
--- a/Carpet/CarpetLib/interface.ccl
+++ b/Carpet/CarpetLib/interface.ccl
@@ -38,6 +38,9 @@ uses include header: CarpetTimers.hh
 uses include header: loopcontrol.h
 uses include header: vectors.h
 
+# DGFE
+uses include header: hrscc.hh
+
 
 
 # Return a pointer to an unmodifiable C string
@@ -69,4 +72,3 @@ CCTK_INT FUNCTION GetBoundarySpecification \
    CCTK_INT OUT ARRAY is_staggered, \
    CCTK_INT OUT ARRAY shiftout)
 USES FUNCTION GetBoundarySpecification
-
diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl
index b4421f097..2ad563960 100644
--- a/Carpet/CarpetLib/param.ccl
+++ b/Carpet/CarpetLib/param.ccl
@@ -27,6 +27,12 @@ BOOLEAN commstate_verbose "Print debug info from the commstate class" STEERABLE=
 
 
 
+BOOLEAN use_dgfe "Use DGFE operators instead of Lagrange operators"
+{
+} "no"
+
+
+
 BOOLEAN output_bboxes "Output bounding box information to the screen" STEERABLE=always
 {
 } "no"
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 24534b376..1ee6c94e7 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -711,6 +711,8 @@ transfer_prolongate (data const * const src,
                      ibbox const & srcbox,
                      int const order_space)
 {
+  DECLARE_CCTK_PARAMETERS;
+  
   static Timer total ("prolongate");
   total.start ();
   
@@ -773,6 +775,17 @@ transfer_prolongate (data const * const src,
       break;
     }
     case cell_centered: {
+      if (use_dgfe) {
+        call_operator<T>(prolongate_3d_dgfe_rf2<T,5>,
+                         static_cast<T const *>(src->storage()),
+                         src->shape(),
+                         static_cast<T *>(this->storage()),
+                         this->shape(),
+                         src->extent(),
+                         this->extent(),
+                         srcbox, dstbox, NULL);
+        break;
+      }
       static
         void (* the_operators[]) (T const * restrict const src,
                                   ivect3 const & restrict srcext,
@@ -1084,6 +1097,8 @@ transfer_restrict (data const * const src,
                    islab const * restrict const slabinfo,
                    int const /*order_space*/)
 {
+  DECLARE_CCTK_PARAMETERS;
+  
   static Timer total ("restrict");
   total.start ();
   
@@ -1119,6 +1134,17 @@ transfer_restrict (data const * const src,
       ibbox const& dstbox = this->extent();
       
       if (all(is_centered == ivect(1,1,1))) {
+        if (use_dgfe) {
+          call_operator<T>(restrict_3d_dgfe_rf2<T,5>,
+                           static_cast<T const *>(src->storage()),
+                           src->shape(),
+                           static_cast<T *>(this->storage()),
+                           this->shape(),
+                           srcbox,
+                           dstbox,
+                           srcregbox, dstregbox, NULL);
+          break;
+        }
         call_operator<T> (& restrict_3d_cc_rf2,
                           static_cast <T const *> (src->storage()),
                           src->shape(),
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index 4f5c36a47..84468a23c 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -35,16 +35,18 @@ SRCS =	balance.cc				\
 	restrict_3d_cc_rf2.cc			\
 	restrict_3d_rf2.cc			\
 	restrict_3d_vc_rf2.cc			\
+	restrict_3d_dgfe_rf2.cc			\
 	restrict_4d_rf2.cc			\
 	prolongate_3d_rf2.cc			\
 	prolongate_3d_cc_rf2.cc			\
+	prolongate_3d_cc_eno_rf2.cc		\
 	prolongate_3d_o5_monotone_rf2.cc	\
 	prolongate_3d_real8_eno.F90		\
 	prolongate_3d_real8_tvd.F90		\
 	prolongate_3d_cc_real8_tvd.F90		\
 	prolongate_3d_real8_weno.F90		\
-	prolongate_4d_o1_rf2.cc			\
-	prolongate_3d_cc_eno_rf2.cc
+	prolongate_3d_dgfe_rf2.cc		\
+	prolongate_4d_o1_rf2.cc
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/Carpet/CarpetLib/src/operator_prototypes_3d.hh b/Carpet/CarpetLib/src/operator_prototypes_3d.hh
index 3cf3cef11..dc775c271 100644
--- a/Carpet/CarpetLib/src/operator_prototypes_3d.hh
+++ b/Carpet/CarpetLib/src/operator_prototypes_3d.hh
@@ -227,6 +227,32 @@ namespace CarpetLib {
                             void * extraargs);
   
   
+  template <typename T, int ORDER>
+  void
+  prolongate_3d_cc_eno_rf2 (T const * restrict const src,
+                            ivect3 const & restrict srcext,
+                            T * restrict const dst,
+                            ivect3 const & restrict dstext,
+                            ibbox3 const & restrict srcbbox,
+                            ibbox3 const & restrict dstbbox,
+                            ibbox3 const & restrict regbbox);
+  
+  
+  
+  template <typename T, int ORDER>
+  void
+  prolongate_3d_dgfe_rf2 (T const * restrict const src,
+                          ivect3 const & restrict srcext,
+                          T * restrict const dst,
+                          ivect3 const & restrict dstext,
+                          ibbox3 const & restrict srcbbox,
+                          ibbox3 const & restrict dstbbox,
+                          ibbox3 const & restrict srcregbbox,
+                          ibbox3 const & restrict dstregbbox,
+                          void * extraargs);
+  
+  
+  
   template <typename T>
   void
   restrict_3d_rf2 (T const * restrict const src,
@@ -365,6 +391,20 @@ namespace CarpetLib {
   
   
   
+  template <typename T, int ORDER>
+  void
+  restrict_3d_dgfe_rf2 (T const * restrict const src,
+                        ivect3 const & restrict srcext,
+                        T * restrict const dst,
+                        ivect3 const & restrict dstext,
+                        ibbox3 const & restrict srcbbox,
+                        ibbox3 const & restrict dstbbox,
+                        ibbox3 const & restrict srcregbbox,
+                        ibbox3 const & restrict dstregbbox,
+                        void * extraargs);
+  
+  
+  
 } // namespace CarpetLib
 
 
diff --git a/Carpet/CarpetLib/src/prolongate_3d_dgfe_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_dgfe_rf2.cc
new file mode 100644
index 000000000..87a219007
--- /dev/null
+++ b/Carpet/CarpetLib/src/prolongate_3d_dgfe_rf2.cc
@@ -0,0 +1,309 @@
+#include <cctk.h>
+#include <cctk_Parameters.h>
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstdlib>
+
+#include <hrscc.hh>
+
+#include "operator_prototypes_3d.hh"
+#include "typeprops.hh"
+
+using namespace std;
+#ifdef HRSCC_HH
+using namespace hrscc;
+#endif
+
+
+
+namespace CarpetLib {
+  
+#define SRCIND3(i,j,k) ptrdiff_t(index3(i, j, k, srciext, srcjext, srckext))
+#define DSTIND3(i,j,k) ptrdiff_t(index3(i, j, k, dstiext, dstjext, dstkext))
+  
+  
+  
+  template<typename T, int ORDER>
+  void
+  prolongate_3d_dgfe_rf2(T const *restrict const src,
+                         ivect3 const& restrict srcext,
+                         T *restrict const dst,
+                         ivect3 const& restrict dstext,
+                         ibbox3 const& restrict srcbbox,
+                         ibbox3 const& restrict dstbbox,
+                         ibbox3 const& restrict,
+                         ibbox3 const& restrict regbbox,
+                         void* const extraargs)
+  {
+    assert(not extraargs);
+    
+    static_assert(ORDER>=0, "ORDER must be non-negative");
+    
+    
+    
+    if (any(srcbbox.stride() <= regbbox.stride() or
+            dstbbox.stride() != regbbox.stride()))
+    {
+      CCTK_WARN(0, "Internal error: strides disagree");
+    }
+    
+    if (any(srcbbox.stride() != reffact2 * dstbbox.stride())) {
+      CCTK_WARN(0, "Internal error: source strides are not twice the destination strides");
+    }
+    
+    if (any(dstbbox.stride() % 2 != 0)) {
+      CCTK_WARN(0, "Internal error: destination strides are not even");
+    }
+    
+    // This could be handled, but is likely to point to an error
+    // elsewhere
+    if (regbbox.empty()) {
+      CCTK_WARN(0, "Internal error: region extent is empty");
+    }
+    
+    
+    
+    bvect3 const is_upper_face =
+      (regbbox.lower() - srcbbox.lower() + regbbox.stride() / 2) %
+      srcbbox.stride() != 0;
+    
+    ivect3 const regext = regbbox.shape() / regbbox.stride();
+    assert(all((regbbox.lower() - srcbbox.lower() +
+                either(is_upper_face, -1, +1) * regbbox.stride() / 2) %
+               srcbbox.stride() == 0));
+    ivect3 const srcoff =
+      (regbbox.lower()- srcbbox.lower() +
+       either(is_upper_face, -1, +1) * regbbox.stride() / 2) /
+      srcbbox.stride();
+    assert(all((regbbox.lower() - dstbbox.lower()) % regbbox.stride() == 0));
+    ivect3 const dstoff =
+      (regbbox.lower() - dstbbox.lower()) / regbbox.stride();
+    
+    
+    
+    if (not regbbox.is_contained_in(srcbbox) or
+        not regbbox.is_contained_in(dstbbox))
+    {
+      cerr << "ORDER=" << ORDER << "\n"
+           << "regbbox=" << regbbox << "\n"
+           << "dstbbox=" << dstbbox << "\n"
+           << "regbbox=" << regbbox << "\n"
+           << "srcbbox=" << srcbbox << "\n";
+      CCTK_WARN(0, "Internal error: region extent is not contained in array extent");
+    }
+    
+    
+    
+    ptrdiff_t const srciext = srcext[0];
+    ptrdiff_t const srcjext = srcext[1];
+    ptrdiff_t const srckext = srcext[2];
+    
+    ptrdiff_t const dstiext = dstext[0];
+    ptrdiff_t const dstjext = dstext[1];
+    ptrdiff_t const dstkext = dstext[2];
+    
+    ptrdiff_t const regiext = regext[0];
+    ptrdiff_t const regjext = regext[1];
+    ptrdiff_t const regkext = regext[2];
+    
+    ptrdiff_t const srcioff = srcoff[0];
+    ptrdiff_t const srcjoff = srcoff[1];
+    ptrdiff_t const srckoff = srcoff[2];
+    
+    ptrdiff_t const dstioff = dstoff[0];
+    ptrdiff_t const dstjoff = dstoff[1];
+    ptrdiff_t const dstkoff = dstoff[2];
+    
+    
+    
+    if (regext[0] == 1) {
+      // 2D prolongation on x face
+      
+      assert(not is_upper_face[1]);
+      assert(not is_upper_face[2]);
+      
+      // Ensure we traverse an even integer number of elements
+      assert(regext[1] % (2*(ORDER+1)) == 0);
+      assert(regext[2] % (2*(ORDER+1)) == 0);
+      
+      int const srcdi = 1;      // 2d face
+      int const srcdj = SRCIND3(0,1,0) - SRCIND3(0,0,0);
+      int const srcdk = SRCIND3(0,0,1) - SRCIND3(0,0,0);
+      int const dstdi = 1;      // 2d face
+      int const dstdj = DSTIND3(0,1,0) - DSTIND3(0,0,0);
+      int const dstdk = DSTIND3(0,0,1) - DSTIND3(0,0,0);
+      int const srcstr2d[2] = {srcdj, srcdk};
+      int const dststr2d[2] = {dstdj, dstdk};
+      
+      // Loop over fine region
+      for (ptrdiff_t k=0; k<regkext; k+=2*(ORDER+1)) {
+        for (ptrdiff_t j=0; j<regjext; j+=2*(ORDER+1)) {
+          ptrdiff_t const i=0;
+#ifdef HRSCC_HH
+          GLLElement<ORDER>::prolongate_2D
+            (&src[SRCIND3(srcioff+i, srcjoff+j, srckoff+k)], srcstr2d,
+             &dst[DSTIND3(dstioff+2*i, dstjoff+2*j, dstkoff+2*k)], dststr2d);
+#else
+          // HRSCCore is not available
+          assert(0);
+#endif
+        }
+      }
+      
+    } else if (regext[1] == 1) {
+      // 2D prolongation on y face
+      
+      assert(not is_upper_face[0]);
+      assert(not is_upper_face[2]);
+      
+      // Ensure we traverse an even integer number of elements
+      assert(regext[0] % (2*(ORDER+1)) == 0);
+      assert(regext[2] % (2*(ORDER+1)) == 0);
+      
+      // int const srcdi = SRCIND3(1,0,0) - SRCIND3(0,0,0);
+      int const srcdi = 1;
+      assert(srcdi == SRCIND3(1,0,0) - SRCIND3(0,0,0));
+      int const srcdj = 1;      // 2d face
+      int const srcdk = SRCIND3(0,0,1) - SRCIND3(0,0,0);
+      // int const dstdi = DSTIND3(1,0,0) - DSTIND3(0,0,0);
+      int const dstdi = 1;
+      assert(dstdi == DSTIND3(1,0,0) - DSTIND3(0,0,0));
+      int const dstdj = 1;      // 2d face
+      int const dstdk = DSTIND3(0,0,1) - DSTIND3(0,0,0);
+      int const srcstr2d[2]= {srcdi, srcdk};
+      int const dststr2d[2]= {dstdi, dstdk};
+      
+      // Loop over fine region
+      for (ptrdiff_t k=0; k<regkext; k+=2*(ORDER+1)) {
+        ptrdiff_t const j=0;
+        for (ptrdiff_t i=0; i<regiext; i+=2*(ORDER+1)) {
+#ifdef HRSCC_HH
+          GLLElement<ORDER>::prolongate_2D
+            (&src[SRCIND3(srcioff+i, srcjoff+j, srckoff+k)], srcstr2d,
+             &dst[DSTIND3(dstioff+2*i, dstjoff+2*j, dstkoff+2*k)], dststr2d);
+#else
+          // HRSCCore is not available
+          assert(0);
+#endif
+        }
+      }
+      
+    } else if (regext[2] == 1) {
+      // 2D prolongation on z face
+      
+      assert(not is_upper_face[0]);
+      assert(not is_upper_face[1]);
+      
+      // Ensure we traverse an even integer number of elements
+      assert(regext[0] % (2*(ORDER+1)) == 0);
+      assert(regext[1] % (2*(ORDER+1)) == 0);
+      
+      // int const srcdi = SRCIND3(1,0,0) - SRCIND3(0,0,0);
+      int const srcdi = 1;
+      assert(srcdi == SRCIND3(1,0,0) - SRCIND3(0,0,0));
+      int const srcdj = SRCIND3(0,1,0) - SRCIND3(0,0,0);
+      int const srcdk = 1;      // 2d face
+      // int const dstdi = DSTIND3(1,0,0) - DSTIND3(0,0,0);
+      int const dstdi = 1;
+      assert(dstdi == DSTIND3(1,0,0) - DSTIND3(0,0,0));
+      int const dstdj = DSTIND3(0,1,0) - DSTIND3(0,0,0);
+      int const dstdk = 1;      // 2d face
+      int const srcstr2d[2]= {srcdi, srcdj};
+      int const dststr2d[2]= {dstdi, dstdj};
+      
+      // Loop over fine region
+      ptrdiff_t const k=0;
+      for (ptrdiff_t j=0; j<regjext; j+=2*(ORDER+1)) {
+        for (ptrdiff_t i=0; i<regiext; i+=2*(ORDER+1)) {
+#ifdef HRSCC_HH
+          GLLElement<ORDER>::prolongate_2D
+            (&src[SRCIND3(srcioff+i, srcjoff+j, srckoff+k)], srcstr2d,
+             &dst[DSTIND3(dstioff+2*i, dstjoff+2*j, dstkoff+2*k)], dststr2d);
+#else
+          // HRSCCore is not available
+          assert(0);
+#endif
+        }
+      }
+      
+    } else {
+      // 3D prolongation
+      
+      assert(not any(is_upper_face));
+      
+      // Ensure we traverse an even integer number of elements
+      assert(all(regext % (2*(ORDER+1)) == 0));
+      
+      // int const srcdi = SRCIND3(1,0,0) - SRCIND3(0,0,0);
+      int const srcdi = 1;
+      assert(srcdi == SRCIND3(1,0,0) - SRCIND3(0,0,0));
+      int const srcdj = SRCIND3(0,1,0) - SRCIND3(0,0,0);
+      int const srcdk = SRCIND3(0,0,1) - SRCIND3(0,0,0);
+      // int const dstdi = DSTIND3(1,0,0) - DSTIND3(0,0,0);
+      int const dstdi = 1;
+      assert(dstdi == DSTIND3(1,0,0) - DSTIND3(0,0,0));
+      int const dstdj = DSTIND3(0,1,0) - DSTIND3(0,0,0);
+      int const dstdk = DSTIND3(0,0,1) - DSTIND3(0,0,0);
+      int const srcstr[3] = {srcdi, srcdj, srcdk};
+      int const dststr[3] = {dstdi, dstdj, dstdk};
+      
+      // Loop over fine region
+      for (ptrdiff_t k=0; k<regkext; k+=2*(ORDER+1)) {
+        for (ptrdiff_t j=0; j<regjext; j+=2*(ORDER+1)) {
+          for (ptrdiff_t i=0; i<regiext; i+=2*(ORDER+1)) {
+#ifdef HRSCC_HH
+            GLLElement<ORDER>::prolongate_full
+              (&src[SRCIND3(srcioff+i, srcjoff+j, srckoff+k)], srcstr,
+               &dst[DSTIND3(dstioff+2*i, dstjoff+2*j, dstkoff+2*k)], dststr);
+#else
+            // HRSCCore is not available
+            assert(0);
+#endif
+          }
+        }
+      }
+      
+    }
+  }
+  
+  
+  
+#define TYPECASE(N,T)                                           \
+                                                                \
+  template                                                      \
+  void                                                          \
+  prolongate_3d_dgfe_rf2<T,5>(T const *restrict const src,      \
+                              ivect3 const& restrict srcext,    \
+                              T *restrict const dst,            \
+                              ivect3 const& restrict dstext,    \
+                              ibbox3 const& restrict srcbbox,   \
+                              ibbox3 const& restrict dstbbox,   \
+                              ibbox3 const& restrict,           \
+                              ibbox3 const& restrict regbbox,   \
+                              void *extraargs);
+  
+#define CARPET_NO_INT
+#define CARPET_NO_COMPLEX
+#include "typecase.hh"
+#undef TYPECASE
+  
+  
+  
+  template<>
+  void
+  prolongate_3d_dgfe_rf2<CCTK_COMPLEX,5>(CCTK_COMPLEX const *restrict const src,
+                                         ivect3 const& restrict srcext,
+                                         CCTK_COMPLEX *restrict const dst,
+                                         ivect3 const& restrict dstext,
+                                         ibbox3 const& restrict srcbbox,
+                                         ibbox3 const& restrict dstbbox,
+                                         ibbox3 const& restrict,
+                                         ibbox3 const& restrict regbbox,
+                                         void *extraargs)
+  {
+    assert(0);
+  }
+  
+} // namespace CarpetLib
diff --git a/Carpet/CarpetLib/src/restrict_3d_dgfe_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_dgfe_rf2.cc
new file mode 100644
index 000000000..42a087692
--- /dev/null
+++ b/Carpet/CarpetLib/src/restrict_3d_dgfe_rf2.cc
@@ -0,0 +1,180 @@
+#include <cctk.h>
+#include <cctk_Parameters.h>
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstdlib>
+
+#include <hrscc.hh>
+
+#include "operator_prototypes_3d.hh"
+#include "typeprops.hh"
+
+using namespace std;
+#ifdef HRSCC_HH
+using namespace hrscc;
+#endif
+
+
+
+namespace CarpetLib {
+  
+#define SRCIND3(i,j,k) ptrdiff_t(index3(i, j, k, srciext, srcjext, srckext))
+#define DSTIND3(i,j,k) ptrdiff_t(index3(i, j, k, dstiext, dstjext, dstkext))
+  
+  
+  
+  template<typename T, int ORDER>
+  void
+  restrict_3d_dgfe_rf2(T const *restrict const src,
+                       ivect3 const& restrict srcext,
+                       T *restrict const dst,
+                       ivect3 const& restrict dstext,
+                       ibbox3 const& restrict srcbbox,
+                       ibbox3 const& restrict dstbbox,
+                       ibbox3 const& restrict,
+                       ibbox3 const& restrict regbbox,
+                       void *const extraargs)
+  {
+    assert(not extraargs);
+    
+    static_assert(ORDER>=0, "ORDER must be non-negative");
+    
+    
+    
+    if (any(srcbbox.stride() >= regbbox.stride() or
+            dstbbox.stride() != regbbox.stride()))
+    {
+      CCTK_WARN(0, "Internal error: strides disagree");
+    }
+    
+    if (any(reffact2 * srcbbox.stride() != dstbbox.stride())) {
+      CCTK_WARN(0, "Internal error: destination strides are not twice the source strides");
+    }
+    
+    // This could be handled, but is likely to point to an error
+    // elsewhere
+    if (regbbox.empty()) {
+      CCTK_WARN(0, "Internal error: region extent is empty");
+    }
+    
+    if (not regbbox.expanded_for(srcbbox).is_contained_in(srcbbox) or
+        not regbbox.is_contained_in(dstbbox))
+    {
+      CCTK_WARN(0, "Internal error: region extent is not contained in array extent");
+    }
+    
+    
+    
+    ivect3 const regext = regbbox.shape() / regbbox.stride();
+    assert(all(srcbbox.stride() % 2 == 0));
+    assert(all((regbbox.lower() - srcbbox.lower() - srcbbox.stride() / 2) %
+               srcbbox.stride() == 0));
+    ivect3 const srcoff =
+      (regbbox.lower() - srcbbox.lower() - srcbbox.stride() / 2) /
+      srcbbox.stride();
+    assert(all((regbbox.lower() - dstbbox.lower()) % dstbbox.stride() == 0));
+    ivect3 const dstoff =
+      (regbbox.lower() - dstbbox.lower()) / dstbbox.stride();
+    
+    
+    
+    ptrdiff_t const srciext = srcext[0];
+    ptrdiff_t const srcjext = srcext[1];
+    ptrdiff_t const srckext = srcext[2];
+    
+    ptrdiff_t const dstiext = dstext[0];
+    ptrdiff_t const dstjext = dstext[1];
+    ptrdiff_t const dstkext = dstext[2];
+    
+    ptrdiff_t const regiext = regext[0];
+    ptrdiff_t const regjext = regext[1];
+    ptrdiff_t const regkext = regext[2];
+    
+    ptrdiff_t const srcioff = srcoff[0];
+    ptrdiff_t const srcjoff = srcoff[1];
+    ptrdiff_t const srckoff = srcoff[2];
+    
+    ptrdiff_t const dstioff = dstoff[0];
+    ptrdiff_t const dstjoff = dstoff[1];
+    ptrdiff_t const dstkoff = dstoff[2];
+    
+    
+    
+    // int const srcdi = SRCIND3(1,0,0) - SRCIND3(0,0,0);
+    int const srcdi = 1;
+    assert(srcdi == SRCIND3(1,0,0) - SRCIND3(0,0,0));
+    int const srcdj = SRCIND3(0,1,0) - SRCIND3(0,0,0);
+    int const srcdk = SRCIND3(0,0,1) - SRCIND3(0,0,0);
+    
+    // int const dstdi = DSTIND3(1,0,0) - DSTIND3(0,0,0);
+    int const dstdi = 1;
+    assert(dstdi == DSTIND3(1,0,0) - DSTIND3(0,0,0));
+    int const dstdj = DSTIND3(0,1,0) - DSTIND3(0,0,0);
+    int const dstdk = DSTIND3(0,0,1) - DSTIND3(0,0,0);
+    
+    int const srcstr[3] = {srcdi, srcdj, srcdk};
+    int const dststr[3] = {dstdi, dstdj, dstdk};
+    
+    
+    
+    // Ensure we traverse an integer number of elements
+    assert(all(regext % (ORDER+1) == 0));
+    
+    // Loop over coarse region
+    for (ptrdiff_t k=0; k<regkext; k+=ORDER+1) {
+      for (ptrdiff_t j=0; j<regjext; j+=ORDER+1) {
+        for (ptrdiff_t i=0; i<regiext; i+=ORDER+1) {
+#ifdef HRSCC_HH
+          GLLElement<ORDER>::restrict_full
+            (&src[SRCIND3(srcioff+2*i, srcjoff+2*j, srckoff+2*k)], srcstr,
+             &dst[DSTIND3(dstioff+i, dstjoff+j, dstkoff+k)], dststr);
+#else
+          // HRSCCore is not available
+          assert(0);
+#endif
+        }
+      }
+    }
+    
+  }
+  
+  
+  
+#define TYPECASE(N,T)                                          \
+  template                                                     \
+  void                                                         \
+  restrict_3d_dgfe_rf2<T,5>(T const *restrict const src,       \
+                            ivect3 const& restrict srcext,     \
+                            T *restrict const dst,             \
+                            ivect3 const& restrict dstext,     \
+                            ibbox3 const& restrict srcbbox,    \
+                            ibbox3 const& restrict dstbbox,    \
+                            ibbox3 const& restrict,            \
+                            ibbox3 const& restrict regbbox,    \
+                            void *extraargs);
+  
+#define CARPET_NO_INT
+#define CARPET_NO_COMPLEX
+#include "typecase.hh"
+#undef TYPECASE
+  
+  
+  
+  template<>
+  void
+  restrict_3d_dgfe_rf2<CCTK_COMPLEX,5>(CCTK_COMPLEX const *restrict const src,
+                                       ivect3 const& restrict srcext,
+                                       CCTK_COMPLEX *restrict const dst,
+                                       ivect3 const& restrict dstext,
+                                       ibbox3 const& restrict srcbbox,
+                                       ibbox3 const& restrict dstbbox,
+                                       ibbox3 const& restrict,
+                                       ibbox3 const& restrict regbbox,
+                                       void *extraargs)
+  {
+    assert(0);
+  }
+  
+} // namespace CarpetLib