First cut at the interpolator. Tentatively works on unigrid and for

First cut at the interpolator.  Tentatively works on unigrid and for
singlepatch mode.

darcs-hash:20030430103731-07bb3-dad06ae817a997cfd4f1beefc4065ca1ecc12342.gz

4 files changed, 75 insertions, 641 deletions

diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc
index 8b48e9c54..47d5b9cb8 100644
--- a/Carpet/CarpetInterp/src/interp.cc
+++ b/Carpet/CarpetInterp/src/interp.cc
@@ -1,30 +1,22 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.36 2004/09/13 14:36:13 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.1 2003/04/30 12:37:31 schnetter Exp $
 
 #include <assert.h>
-#include <math.h>
 
-#include <algorithm>
+#include <complex>
 #include <vector>
 
 #include <mpi.h>
 
 #include "cctk.h"
 
-#include "util_ErrorCodes.h"
-#include "util_Table.h"
+#include "Carpet/CarpetLib/src/vect.hh"
 
-#include "bbox.hh"
-#include "data.hh"
-#include "defs.hh"
-#include "ggf.hh"
-#include "vect.hh"
-
-#include "carpet.hh"
+#include "Carpet/Carpet/src/carpet.hh"
 
 #include "interp.hh"
 
 extern "C" {
-  static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.36 2004/09/13 14:36:13 schnetter Exp $";
+  static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.1 2003/04/30 12:37:31 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetInterp_interp_cc);
 }
 
@@ -39,84 +31,6 @@ namespace CarpetInterp {
   typedef vect<int,dim> ivect;
   typedef vect<CCTK_REAL,dim> rvect;
   
-  typedef bbox<int,dim> ibbox;
-  
-  
-  
-#define ind_rc(m,rl,c) ind_rc_(m,rl,minrl,maxrl,c,maxncomps,vhh)
-  static inline int ind_rc_(int const m,
-                            int const rl, int const minrl, int const maxrl,
-                            int const c, int const maxncomps,
-                            vector<gh<dim> *> const hh)
-  {
-    assert (m>=0 && m<maps);
-    assert (rl>=minrl && rl<maxrl);
-    assert (minrl>=0 && maxrl<=hh.at(m)->reflevels());
-    assert (c>=0 && c<maxncomps);
-    assert (maxncomps>=0 && maxncomps<=hh.at(m)->components(rl));
-    int const ind = (rl-minrl) * maxncomps + c;
-    assert (ind>=0 && ind < (maxrl-minrl) * maxncomps);
-    return ind;
-  }
-  
-#define ind_prc(p,m,rl,c) \
-  ind_prc_(p,nprocs,m,rl,minrl,maxrl,c,maxncomps,cgh,vhh)
-  static inline int ind_prc_(int const p, int const nprocs,
-                             int const m,
-                             int const rl, int const minrl, int const maxrl,
-                             int const c, int const maxncomps,
-                             cGH const * const cgh,
-                             vector<gh<dim> *> const hh)
-  {
-    assert (p>=0 && p<nprocs);
-    assert (nprocs==CCTK_nProcs(cgh));
-    assert (m>=0 && m<maps);
-    assert (rl>=minrl && rl<maxrl);
-    assert (minrl>=0 && maxrl<=hh.at(m)->reflevels());
-    assert (c>=0 && c<maxncomps);
-    assert (maxncomps>=0 && maxncomps<=hh.at(m)->components(rl));
-    int const ind = (p * (maxrl-minrl) + (rl-minrl)) * maxncomps + c;
-    assert (ind>=0 && ind < nprocs * (maxrl-minrl) * maxncomps);
-    return ind;
-  }
-  
-
-  static int GetInterpNumTimelevels(const cGH * const cgh,
-                                    const int group)
-  {
-    assert (group>=0 && group<CCTK_NumGroups());
-    
-    int ierr;
-    
-    cGroup gp;
-    ierr = CCTK_GroupData (group, &gp);
-    assert (!ierr);
-    
-    int interp_ntl = -1;
-
-    const int length = Util_TableGetInt
-      (gp.tagstable, &interp_ntl, "InterpNumTimelevels");
-    if (length == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
-      interp_ntl = -1;
-    } else if (length >= 0) {
-      // all is well - just check they're not asking too much
-      if (interp_ntl > CCTK_ActiveTimeLevelsGI(cgh, group)) {
-        char * const groupname = CCTK_GroupName (group);
-        CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
-                   "The tags table entry \"InterpNumTimelevels\" "
-                   "for group \"%s\" says that %d timelevels should "
-                   "be used for time interpolation.\n However, only %d "
-                   "are active.", 
-                   groupname, interp_ntl, 
-                   CCTK_ActiveTimeLevelsGI(cgh, group));
-      }             
-    } else {
-      assert(0);
-    }
-
-    return interp_ntl;
-    
-  }
   
   
   void CarpetInterpStartup ()
@@ -126,7 +40,7 @@ namespace CarpetInterp {
   
   
   
-  int InterpGridArrays (cGH const * const cgh,
+  int InterpGridArrays (cGH const * const cGH,
                         int const N_dims,
                         int const local_interp_handle,
                         int const param_table_handle,
@@ -140,556 +54,77 @@ namespace CarpetInterp {
                         CCTK_INT const output_array_type_codes [],
                         void * const output_arrays [])
   {
-    if (CCTK_IsFunctionAliased ("SymmetryInterpolate")) {
-      return SymmetryInterpolate
-        (cgh, N_dims,
-         local_interp_handle, param_table_handle, coord_system_handle,
-         N_interp_points, interp_coords_type_code, interp_coords,
-         N_input_arrays, input_array_variable_indices,
-         N_output_arrays, output_array_type_codes, output_arrays);
-    } else {
-      return Carpet_DriverInterpolate
-        (cgh, N_dims,
-         local_interp_handle, param_table_handle, coord_system_handle,
-         N_interp_points, interp_coords_type_code, interp_coords,
-         N_input_arrays, input_array_variable_indices,
-         N_output_arrays, output_array_type_codes, output_arrays);
-    }
-  }
-  
-  
-  
-  extern "C" CCTK_INT
-  Carpet_DriverInterpolate (CCTK_POINTER_TO_CONST const cgh_,
-                            CCTK_INT const N_dims,
-                            CCTK_INT const local_interp_handle,
-                            CCTK_INT const param_table_handle,
-                            CCTK_INT const coord_system_handle,
-                            CCTK_INT const N_interp_points,
-                            CCTK_INT const interp_coords_type_code,
-                            CCTK_POINTER_TO_CONST const interp_coords [],
-                            CCTK_INT const N_input_arrays,
-                            CCTK_INT const input_array_variable_indices [],
-                            CCTK_INT const N_output_arrays,
-                            CCTK_INT const output_array_type_codes [],
-                            CCTK_POINTER const output_arrays [])
-  {
-    cGH const * const cgh = static_cast<cGH const *> (cgh_);
     int ierr;
     
-    assert (cgh);
+    assert (cGH);
     assert (N_dims == dim);
     
-    
-    
-    if (is_meta_mode()) {
-      CCTK_WARN (0, "It is not possible to interpolate in meta mode");
-    }
-    
-    bool const want_global_mode = is_global_mode();
-    
-    
-    
-    // Get some information
-    MPI_Comm const comm = CarpetMPIComm ();
-    int const myproc = CCTK_MyProc (cgh);
-    int const nprocs = CCTK_nProcs (cgh);
-    assert (myproc>=0 && myproc<nprocs);
-    
-    // Multiple maps are not supported
-    // (because we don't know how to select a map)
-    assert (maps == 1);
-    const int m = 0;
-    
-    int const minrl = want_global_mode ? 0                      : reflevel;
-    int const maxrl = want_global_mode ? vhh.at(m)->reflevels() : reflevel+1;
-    
-    // Multiple convergence levels are not supported
-    assert (mglevels == 1);
-    int const ml = 0;
-    
-    int maxncomps = 0;
-    for (int rl=minrl; rl<maxrl; ++rl) {
-      maxncomps = max(maxncomps, vhh.at(m)->components(rl));
-    }
-    
-    
-    
-    // Find the time interpolation order
-    int partype;
-    void const * const parptr
-      = CCTK_ParameterGet ("prolongation_order_time", "Carpet", &partype);
-    assert (parptr);
-    assert (partype == PARAMETER_INTEGER);
-    int const prolongation_order_time = * (CCTK_INT const *) parptr;
-    
-    CCTK_REAL const current_time = cgh->cctk_time / cgh->cctk_delta_time;
-    
-    
-    
-    // Find out about the coordinates: origin and delta for the Carpet
-    // grid indices
-#if 0
     const char * coord_system_name
       = CCTK_CoordSystemName (coord_system_handle);
     assert (coord_system_name);
-    rvect lower, upper, delta;
+    rvect lower, upper;
     for (int d=0; d<dim; ++d) {
       ierr = CCTK_CoordRange
-        (cgh, &lower[d], &upper[d], d+1, 0, coord_system_name);
+        (cGH, &lower[d], &upper[d], d+1, 0, coord_system_name);
       assert (!ierr);
-      const ibbox & baseext = vhh.at(m)->baseextent;
-      delta[d]
-        = (upper[d] - lower[d]) / (baseext.upper()[d] - baseext.lower()[d]);
-    }
-#else
-    rvect const lower = rvect::ref(cgh->cctk_origin_space);
-    rvect const delta = rvect::ref(cgh->cctk_delta_space) / maxreflevelfact;
-    rvect const upper = lower + delta * (vhh.at(m)->baseextent.upper() - vhh.at(m)->baseextent.lower());
-#endif
-    
-    assert (N_interp_points >= 0);
-    assert (interp_coords);
-    for (int d=0; d<dim; ++d) {
-      assert (N_interp_points==0 || interp_coords[d]);
     }
     
-    assert (N_output_arrays >= 0);
-    if (N_interp_points > 0) {
-      assert (output_arrays);
-      for (int j=0; j<N_output_arrays; ++j) {
-        assert (output_arrays[j]);
-      }
-    }
-    
-    
+    const ivect lbnd (cGH->cctk_lbnd);
+    const ivect lsh  (cGH->cctk_lsh );
+    const ivect gsh  (cGH->cctk_gsh );
     
-    // Assign interpolation points to components
-    vector<int> rlev (N_interp_points); // refinement level of point n
-    vector<int> home (N_interp_points); // component of point n
-    vector<int> homecnts ((maxrl-minrl) * maxncomps); // number of points in component c
-    
-    for (int n=0; n<N_interp_points; ++n) {
-      
-      // Find the grid point closest to the interpolation point
-      assert (interp_coords_type_code == CCTK_VARIABLE_REAL);
-      rvect pos;
-      for (int d=0; d<dim; ++d) {
-        pos[d] = static_cast<CCTK_REAL const *>(interp_coords[d])[n];
-      }
-      
-      // Find the component that this grid point belongs to
-      rlev.at(n) = -1;
-      home.at(n) = -1;
-      if (all(pos>=lower && pos<=upper)) {
-        for (int rl=maxrl-1; rl>=minrl; --rl) {
-          
-          const int fact = maxreflevelfact / ipow(reffact, rl) * ipow(mgfact, mglevel);
-          ivect const ipos = ivect(floor((pos - lower) / (delta * fact) + 0.5)) * fact;
-          assert (all(ipos % vhh.at(m)->bases.at(rl).at(ml).stride() == 0));
-          
-          // TODO: use something faster than a linear search
-          for (int c=0; c<vhh.at(m)->components(rl); ++c) {
-            if (vhh.at(m)->extents.at(rl).at(c).at(ml).contains(ipos)) {
-              rlev.at(n) = rl;
-              home.at(n) = c;
-              goto found;
-            }
-          }
-        }
-      }
-      CCTK_VWarn (3, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Interpolation point #%d at [%g,%g,%g] is not on any grid patch",
-                  n, pos[0], pos[1], pos[2]);
-      // Map the point (arbitrarily) to the first component of the
-      // coarsest grid
-      rlev.at(n) = minrl;
-      home.at(n) = 0;
-    found:
-      assert (rlev.at(n)>=minrl && rlev.at(n)<maxrl);
-      assert (home.at(n)>=0 && home.at(n)<vhh.at(m)->components(rlev.at(n)));
-      ++ homecnts.at(ind_rc(m, rlev.at(n), home.at(n)));
-      
-    } // for n
-    
-    // Communicate counts
-    vector<int> allhomecnts(nprocs * (maxrl-minrl) * maxncomps);
-    MPI_Allgather
-      (&homecnts   .at(0), (maxrl-minrl) * maxncomps, MPI_INT,
-       &allhomecnts.at(0), (maxrl-minrl) * maxncomps, MPI_INT, comm);
-    
-    
-    
-    // Create coordinate patches
-    vector<data<CCTK_REAL,dim> > allcoords
-      (nprocs * (maxrl-minrl) * maxncomps);
-    for (int p=0; p<nprocs; ++p) {
-      for (int rl=minrl; rl<maxrl; ++rl) {
-        for (int c=0; c<vhh.at(m)->components(rl); ++c) {
-          ivect lo (0);
-          ivect up (1);
-          up[0] = allhomecnts.at(ind_prc(p,m,rl,c));
-          up[1] = dim;
-          ivect str (1);
-          ibbox extent (lo, up-str, str);
-          allcoords.at(ind_prc(p,m,rl,c)).allocate (extent, p);
-        }
-      }
-    }
-    
-    // Fill in local coordinate patches
-    {
-      vector<int> tmpcnts ((maxrl-minrl) * maxncomps);
-      for (int n=0; n<N_interp_points; ++n) {
-        int const rl = rlev.at(n);
-        int const c = home.at(n);
-        assert (rl>=minrl && rl<maxrl);
-        assert (c>=0 && c<vhh.at(m)->components(rl));
-        assert (tmpcnts.at(ind_rc(m,rl,c)) >= 0);
-        assert (tmpcnts.at(ind_rc(m,rl,c)) < homecnts.at(ind_rc(m,rl,c)));
-        assert (dim==3);
-        for (int d=0; d<dim; ++d) {
-          allcoords.at(ind_prc(myproc,m,rl,c))[ivect(tmpcnts.at(ind_rc(m,rl,c)),d,0)]
-            = static_cast<CCTK_REAL const *>(interp_coords[d])[n];
-        }
-        ++ tmpcnts.at(c + (rl-minrl)*maxncomps);
-      }
-      for (int rl=minrl; rl<maxrl; ++rl) {
-        for (int c=0; c<vhh.at(m)->components(rl); ++c) {
-          assert (tmpcnts.at(ind_rc(m,rl,c)) == homecnts.at(ind_rc(m,rl,c)));
-        }
-      }
-    }
-    
-    // Transfer coordinate patches
-    for (comm_state<dim> state; !state.done(); state.step()) {
-      for (int p=0; p<nprocs; ++p) {
-        for (int rl=minrl; rl<maxrl; ++rl) {
-          for (int c=0; c<vhh.at(m)->components(rl); ++c) {
-            allcoords.at(ind_prc(p,m,rl,c)).change_processor
-              (state, vhh.at(m)->processors.at(rl).at(c));
-          }
-        }
-      }
-    }
-    
-    
-    
-    // Create output patches
-    vector<data<CCTK_REAL,dim> > alloutputs
-      (nprocs * (maxrl-minrl) * maxncomps, -1);
-    for (int p=0; p<nprocs; ++p) {
-      for (int rl=minrl; rl<maxrl; ++rl) {
-        for (int c=0; c<vhh.at(m)->components(rl); ++c) {
-          ivect lo (0);
-          ivect up (1);
-          up[0] = allhomecnts.at(ind_prc(p,m,rl,c));
-          up[1] = N_output_arrays;
-          ivect str (1);
-          ibbox extent (lo, up-str, str);
-          alloutputs.at(ind_prc(p,m,rl,c)).allocate
-            (extent, vhh.at(m)->processors.at(rl).at(c));
-        }
-      }
-    }
+    const rvect coord_delta
+      = (upper - lower) / rvect((gsh - 1) * reflevelfact);
+    const rvect coord_origin = lower + rvect(lbnd) * coord_delta;
     
-    
-    
-    //
-    // Do the local interpolation
-    //
-    int overall_ierr = 0;
-    BEGIN_GLOBAL_MODE(cgh) {
-      for (int rl=minrl; rl<maxrl; ++rl) {
-        enter_level_mode (const_cast<cGH*>(cgh), rl);
+    vector<CCTK_INT> input_array_type_codes(N_input_arrays);
+    vector<const void *> input_arrays(N_input_arrays);
+    for (int n=0; n<N_input_arrays; ++n) {
+      if (input_array_variable_indices[n] == -1) {
         
-        // Number of necessary time levels
-        CCTK_REAL const level_time = cgh->cctk_time / cgh->cctk_delta_time;
-        bool const need_time_interp
-          = (fabs(current_time - level_time)
-             > 1e-12 * (fabs(level_time) + fabs(current_time)
-                        + fabs(cgh->cctk_delta_time)));
-        assert (! (! want_global_mode && need_time_interp));
-        int const num_tl = need_time_interp ? prolongation_order_time + 1 : 1;
+        // Ignore this entry
+        input_array_type_codes[n] = -1;
+        input_arrays[n] = 0;
+        
+      } else {
+        
+        int const vi = input_array_variable_indices[n];
+        assert (vi>=0 && vi<CCTK_NumVars());
+        
+        int const gi = CCTK_GroupIndexFromVarI (vi);
+        assert (gi>=0 && gi<CCTK_NumGroups());
+        
+        cGroup group;
+        ierr = CCTK_GroupData (gi, &group);
+        assert (!ierr);
+        assert (group.grouptype == CCTK_GF);
+        assert (group.dim == dim);
+        assert (group.disttype == CCTK_DISTRIB_DEFAULT);
+        assert (group.stagtype == 0); // not staggered
+        
+        input_array_type_codes[n] = group.vartype;
+        input_arrays[n] = CCTK_VarDataPtrI (cGH, 0, vi);
         
-        BEGIN_MAP_LOOP(cgh, CCTK_GF) {
-          BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) {
-            
-            // Find out about the local geometry
-            ivect lsh;
-            rvect coord_origin, coord_delta;
-            for (int d=0; d<dim; ++d) {
-              lsh[d] = cgh->cctk_lsh[d];
-              coord_delta[d] = cgh->cctk_delta_space[d] / cgh->cctk_levfac[d];
-              coord_origin[d] = cgh->cctk_origin_space[d] + (cgh->cctk_lbnd[d] + 1.0 * cgh->cctk_levoff[d] / cgh->cctk_levoffdenom[d]) * coord_delta[d];
-            }
-            
-            
-            
-            // Find out about the grid functions
-            vector<CCTK_INT> input_array_type_codes(N_input_arrays);
-            vector<const void *> input_arrays(N_input_arrays);
-            for (int n=0; n<N_input_arrays; ++n) {
-              if (input_array_variable_indices[n] == -1) {
-                
-                // Ignore this entry
-                input_array_type_codes.at(n) = -1;
-                
-              } else {
-                
-                int const vi = input_array_variable_indices[n];
-                assert (vi>=0 && vi<CCTK_NumVars());
-                
-                int const gi = CCTK_GroupIndexFromVarI (vi);
-                assert (gi>=0 && gi<CCTK_NumGroups());
-                
-                cGroup group;
-                ierr = CCTK_GroupData (gi, &group);
-                assert (!ierr);
-                assert (group.grouptype == CCTK_GF);
-                assert (group.dim == dim);
-                assert (group.disttype == CCTK_DISTRIB_DEFAULT);
-                assert (group.stagtype == 0); // not staggered
-                
-                int interp_ntls = GetInterpNumTimelevels(cgh, gi);
-                // If -ve then key wasn't set; use all timelevels.
-                if (interp_ntls < 0) interp_ntls = group.numtimelevels;
-
-                // TODO: emit better warning
-                if ((num_tl > group.numtimelevels)&&
-                    (interp_ntls > group.numtimelevels)) {
-                  CCTK_VWarn(0, __LINE__,__FILE__,"CarpetInterp",
-                             "Tried to interpolate variable '%s' "
-                             "in time.\nIt has insufficient timelevels "
-                             "(%d are required).",
-                             CCTK_FullName(vi),
-                             min(num_tl,interp_ntls));
-                }
-                assert (group.numtimelevels >= min(num_tl,interp_ntls));
-                
-                input_array_type_codes.at(n) = group.vartype;
-                
-              }
-            } // for input arrays
-            
-            
-            
-              // Work on the data from all processors
-            for (int p=0; p<nprocs; ++p) {
-              assert (allcoords.at(ind_prc(p,m,reflevel,component)).owns_storage());
-              assert (allhomecnts.at(ind_prc(p,m,reflevel,component)) == allcoords.at(ind_prc(p,m,reflevel,component)).shape()[0]);
-              assert (allhomecnts.at(ind_prc(p,m,reflevel,component)) == alloutputs.at(ind_prc(p,m,reflevel,component)).shape()[0]);
-              
-              int const npoints = allhomecnts.at(ind_prc(p,m,reflevel,component));
-              
-              // Do the processor-local interpolation
-              vector<const void *> tmp_interp_coords (dim);
-              for (int d=0; d<dim; ++d) {
-                tmp_interp_coords.at(d) = &allcoords.at(ind_prc(p,m,reflevel,component))[ivect(0,d,0)];
-              }
-              
-              
-              
-              vector<vector<void *> > tmp_output_arrays (num_tl);
-              
-              for (int tl=0; tl<num_tl; ++tl) {
-                
-                for (int n=0; n<N_input_arrays; ++n) {
-
-                  int const vi = input_array_variable_indices[n];
-                  assert (vi>=0 && vi<CCTK_NumVars());
-                
-                  int const gi = CCTK_GroupIndexFromVarI (vi);
-                  assert (gi>=0 && gi<CCTK_NumGroups());
-                
-                  int interp_ntls = GetInterpNumTimelevels(cgh, gi);
-                  // If -ve then key wasn't set; use all timelevels.
-                  if (interp_ntls < 0) interp_ntls = num_tl;
-
-                  if (input_array_variable_indices[n] == -1) {
-                    input_arrays.at(n) = 0;
-                  } else if (tl > interp_ntls-1) {
-                    // Do the interpolation anyway, just from 
-                    // an earlier timelevel.
-                    int const vi = input_array_variable_indices[n];
-                    assert (vi>=0 && vi<CCTK_NumVars());
-                    input_arrays.at(n) = 
-                      CCTK_VarDataPtrI (cgh, interp_ntls-1, vi);
-                  } else {
-                    int const vi = input_array_variable_indices[n];
-                    assert (vi>=0 && vi<CCTK_NumVars());
-                    input_arrays.at(n) = CCTK_VarDataPtrI (cgh, tl, vi);
-                  }
-                } // for input arrays
-                
-                tmp_output_arrays.at(tl).resize (N_output_arrays);
-                for (int j=0; j<N_output_arrays; ++j) {
-                  assert (output_array_type_codes[j] == CCTK_VARIABLE_REAL);
-                  if (need_time_interp) {
-                    tmp_output_arrays.at(tl).at(j) = new CCTK_REAL [npoints];
-                  } else {
-                    tmp_output_arrays.at(tl).at(j) = &alloutputs.at(ind_prc(p,m,reflevel,component))[ivect(0,j,0)];
-                  }
-                }
-                  
-                ierr = CCTK_InterpLocalUniform
-                  (N_dims, local_interp_handle, param_table_handle,
-                   &coord_origin[0], &coord_delta[0],
-                   npoints,
-                   interp_coords_type_code, &tmp_interp_coords[0],
-                   N_input_arrays, &lsh[0],
-                   &input_array_type_codes[0], &input_arrays[0],
-                   N_output_arrays,
-                   output_array_type_codes, &tmp_output_arrays.at(tl)[0]);
-                if (ierr) {
-                  CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                              "The local interpolator returned the error code %d", ierr);
-                }
-                overall_ierr = min(overall_ierr, ierr);
-                
-              } // for tl
-              
-              // Interpolate in time, if necessary
-              if (need_time_interp) {
-                
-                for (int j=0; j<N_output_arrays; ++j) {
-                  
-                  // Find output variable indices
-                  vector<CCTK_INT> operand_indices (N_output_arrays);
-                  ierr = Util_TableGetIntArray
-                    (param_table_handle, N_output_arrays,
-                     &operand_indices.front(), "operand_indices");
-                  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
-                    assert (N_output_arrays == N_input_arrays);
-                    for (int m=0; m<N_output_arrays; ++m) {
-                      operand_indices.at(m) = m;
-                    }
-                  } else {
-                    assert (ierr == N_output_arrays);
-                  }
-                  
-                  // Find input array for this output array
-                  int const n = operand_indices.at(j);
-                  
-                  int const vi = input_array_variable_indices[n];
-                  assert (vi>=0 && vi<CCTK_NumVars());
-                  
-                  int const gi = CCTK_GroupIndexFromVarI (vi);
-                  assert (gi>=0 && gi<CCTK_NumGroups());
-                  
-                  int interp_ntls = GetInterpNumTimelevels(cgh, gi);
-                  // If -ve then key wasn't set; use all timelevels.
-                  if (interp_ntls < 0) interp_ntls = num_tl;
-                  
-                  // Get interpolation times
-                  CCTK_REAL const time = current_time;
-                  vector<CCTK_REAL> times(interp_ntls);
-                  for (int tl=0; tl<interp_ntls; ++tl) {
-                    times.at(tl) = vtt.at(m)->time (-tl, reflevel, mglevel);
-                  }
-                  
-                  // Calculate interpolation weights
-                  vector<CCTK_REAL> tfacs(interp_ntls);
-                  switch (interp_ntls) {
-                    case 1:
-                      // no interpolation
-                      // We have to assume that any GF with one timelevel
-                      // is constant in time!!!
-                      // assert (fabs((time - times.at(0)) / fabs(time + times.at(0) + cgh->cctk_delta_time)) < 1e-12);
-                      tfacs.at(0) = 1.0;
-                      break;
-                    case 2:
-                      // linear (2-point) interpolation
-                      tfacs.at(0) = (time - times.at(1)) / (times.at(0) - times.at(1));
-                      tfacs.at(1) = (time - times.at(0)) / (times.at(1) - times.at(0));
-                      break;
-                    case 3:
-                      // quadratic (3-point) interpolation
-                      tfacs.at(0) = (time - times.at(1)) * (time - times.at(2)) / ((times.at(0) - times.at(1)) * (times.at(0) - times.at(2)));
-                      tfacs.at(1) = (time - times.at(0)) * (time - times.at(2)) / ((times.at(1) - times.at(0)) * (times.at(1) - times.at(2)));
-                      tfacs.at(2) = (time - times.at(0)) * (time - times.at(1)) / ((times.at(2) - times.at(0)) * (times.at(2) - times.at(1)));
-                      break;
-                    default:
-                      assert (0);
-                  }
-                  
-                  // Interpolate
-                  assert (output_array_type_codes[j] == CCTK_VARIABLE_REAL);
-                  for (int k=0; k<npoints; ++k) {
-                    CCTK_REAL & dest = alloutputs.at(ind_prc(p,m,reflevel,component))[ivect(k,j,0)];
-                    dest = 0;
-                    for (int tl=0; tl<interp_ntls; ++tl) {
-                      CCTK_REAL const src = ((CCTK_REAL const *)tmp_output_arrays.at(tl).at(j))[k];
-                      dest += tfacs[tl] * src;
-                    }
-                  }
-                }
-                  
-                for (int tl=0; tl<num_tl; ++tl) {
-                  for (int j=0; j<N_output_arrays; ++j) {
-                    delete [] (CCTK_REAL *)tmp_output_arrays.at(tl).at(j);
-                  }
-                }
-                  
-              } // if need_time_interp
-                
-            } // for processors
-              
-          } END_LOCAL_COMPONENT_LOOP;
-        } END_MAP_LOOP;
-          
-        leave_level_mode (const_cast<cGH*>(cgh));
-      } // for rl
-      
-    } END_GLOBAL_MODE;
-    
-    
-    
-    // Transfer output patches back
-    for (comm_state<dim> state; !state.done(); state.step()) {
-      for (int p=0; p<nprocs; ++p) {
-        for (int rl=minrl; rl<maxrl; ++rl) {
-          for (int c=0; c<vhh.at(m)->components(rl); ++c) {
-            alloutputs.at(ind_prc(p,m,rl,c)).change_processor (state, p);
-          }
-        }
-      }
-    }
-    
-    
-    
-    // Read out local output patches
-    {
-      vector<int> tmpcnts ((maxrl-minrl) * maxncomps);
-      for (int n=0; n<N_interp_points; ++n) {
-        int const rl = rlev.at(n);
-        int const c = home.at(n);
-        for (int j=0; j<N_output_arrays; ++j) {
-          assert (interp_coords_type_code == CCTK_VARIABLE_REAL);
-          assert (alloutputs.at(ind_prc(myproc,m,rl,c)).owns_storage());
-          assert (output_arrays);
-          assert (output_arrays[j]);
-          static_cast<CCTK_REAL *>(output_arrays[j])[n] = alloutputs.at(ind_prc(myproc,m,rl,c))[ivect(tmpcnts.at(ind_rc(m,rl,c)),j,0)];
-        }
-        ++ tmpcnts.at(ind_rc(m,rl,c));
-      }
-      for (int rl=minrl; rl<maxrl; ++rl) {
-        for (int c=0; c<vhh.at(m)->components(rl); ++c) {
-          assert (tmpcnts.at(ind_rc(m,rl,c)) == homecnts.at(ind_rc(m,rl,c)));
-        }
       }
     }
     
-    
-    
-    int global_overall_ierr;
-    MPI_Allreduce
-      (&overall_ierr, &global_overall_ierr, 1, MPI_INT, MPI_MIN, comm);
-    
-    
-    
-    // Done.
-    return global_overall_ierr;
+    ierr = CCTK_InterpLocalUniform (N_dims,
+                                    local_interp_handle,
+                                    param_table_handle,
+                                    &coord_origin[0],
+                                    &coord_delta[0],
+                                    N_interp_points,
+                                    interp_coords_type_code,
+                                    interp_coords,
+                                    N_input_arrays,
+                                    &lsh[0],
+                                    &input_array_type_codes[0],
+                                    &input_arrays[0],
+                                    N_output_arrays,
+                                    output_array_type_codes,
+                                    output_arrays);
+    
+    return ierr;
   }
   
 } // namespace CarpetInterp
diff --git a/Carpet/CarpetInterp/src/interp.h b/Carpet/CarpetInterp/src/interp.h
index 4600c6922..8b349b0da 100644
--- a/Carpet/CarpetInterp/src/interp.h
+++ b/Carpet/CarpetInterp/src/interp.h
@@ -1,4 +1,4 @@
-/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.h,v 1.2 2003/11/05 16:18:38 schnetter Exp $ */
+/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.h,v 1.1 2003/04/30 12:37:31 schnetter Exp $ */
 
 #ifndef CARPETINTERP_H
 #define CARPETINTERP_H
diff --git a/Carpet/CarpetInterp/src/interp.hh b/Carpet/CarpetInterp/src/interp.hh
index db0d1360a..ff81643f6 100644
--- a/Carpet/CarpetInterp/src/interp.hh
+++ b/Carpet/CarpetInterp/src/interp.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.hh,v 1.2 2004/05/21 18:12:23 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.hh,v 1.1 2003/04/30 12:37:31 schnetter Exp $
 
 #ifndef CARPETINTERP_HH
 #define CARPETINTERP_HH
@@ -9,20 +9,19 @@
 
 namespace CarpetInterp {
   
-  int
-  InterpGridArrays (cGH const * const cGH,
-                    int const N_dims,
-                    int const local_interp_handle,
-                    int const param_table_handle,
-                    int const coord_system_handle,
-                    int const N_interp_points,
-                    int const interp_coords_type_code,
-                    void const * const interp_coords [],
-                    int const N_input_arrays,
-                    CCTK_INT const input_array_variable_indices [],
-                    int const N_output_arrays,
-                    CCTK_INT const output_array_type_codes [],
-                    void * const output_arrays []);
+  int InterpGridArrays (cGH const * const cGH,
+                        int const N_dims,
+                        int const local_interp_handle,
+                        int const param_table_handle,
+                        int const coord_system_handle,
+                        int const N_interp_points,
+                        int const interp_coords_type_code,
+                        void const * const interp_coords [],
+                        int const N_input_arrays,
+                        CCTK_INT const input_array_variable_indices [],
+                        int const N_output_arrays,
+                        CCTK_INT const output_array_type_codes [],
+                        void * const output_arrays []);
   
 } // namespace CarpetInterp
 
diff --git a/Carpet/CarpetInterp/src/make.code.defn b/Carpet/CarpetInterp/src/make.code.defn
index d5f82ae59..8511e0de3 100644
--- a/Carpet/CarpetInterp/src/make.code.defn
+++ b/Carpet/CarpetInterp/src/make.code.defn
@@ -1,8 +1,8 @@
 # Main make.code.defn file for thorn CarpetInterp
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/make.code.defn,v 1.1 2003/04/29 14:02:25 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/make.code.defn,v 1.2 2003/04/30 12:37:31 schnetter Exp $
 
 # Source files in this directory
-SRCS = 
+SRCS = interp.cc
 
 # Subdirectories containing source files
 SUBDIRS =