Add parallel (distributed) arrays, and facilities to interpolate to them.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/InterpToArray/trunk@10 6ca9aeac-0e4f-0410-a746-fe4df63e9d0c

4 files changed, 528 insertions, 48 deletions

diff --git a/interface.ccl b/interface.ccl
index 551200d..2a4d687 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -5,6 +5,11 @@ IMPLEMENTS: InterpToArray
 PUBLIC:
 
 REAL scalars[nscalars] TYPE=scalar
+
 REAL arrays1d[narrays1d] TYPE=array DIM=1 DISTRIB=constant SIZE=array1d_npoints_i
 REAL arrays2d[narrays2d] TYPE=array DIM=2 DISTRIB=constant SIZE=array2d_npoints_i,array2d_npoints_j
 REAL arrays3d[narrays3d] TYPE=array DIM=3 DISTRIB=constant SIZE=array3d_npoints_i,array3d_npoints_j,array3d_npoints_k
+
+REAL parrays1d[nparrays1d] TYPE=array DIM=1 DISTRIB=default SIZE=parray1d_npoints_i
+REAL parrays2d[nparrays2d] TYPE=array DIM=2 DISTRIB=default SIZE=parray2d_npoints_i,parray2d_npoints_j
+REAL parrays3d[nparrays3d] TYPE=array DIM=3 DISTRIB=default SIZE=parray3d_npoints_i,parray3d_npoints_j,parray3d_npoints_k
diff --git a/param.ccl b/param.ccl
index aa045d7..83f682a 100644
--- a/param.ccl
+++ b/param.ccl
@@ -52,7 +52,7 @@ INT narrays1d "Number of 1D grid arrays"
   0:100 :: ""
 } 0
 
-STRING array1d_vars[100] "Names of the grid functions that should be interpolated on a line point" STEERABLE=always
+STRING array1d_vars[100] "Names of the grid functions that should be interpolated on a line" STEERABLE=always
 {
   "^$"                                                 :: "do not interpolate"
   "^[A-Za-z][A-Za-z0-9_]*[:][:][A-Za-z][A-Za-z0-9_]*$" :: "grid function name"
@@ -258,3 +258,217 @@ REAL array3d_dz_k "Spacing" STEERABLE=always
 {
   0.0:* :: ""
 } 0.0
+
+
+
+INT nparrays1d "Number of 1D parallel grid arrays"
+{
+  0:100 :: ""
+} 0
+
+STRING parray1d_vars[100] "Names of the grid functions that should be interpolated on a line" STEERABLE=always
+{
+  "^$"                                                 :: "do not interpolate"
+  "^[A-Za-z][A-Za-z0-9_]*[:][:][A-Za-z][A-Za-z0-9_]*$" :: "grid function name"
+} ""
+
+INT parray1d_spacederivs[100] "Space derivative orders for each grid function"
+{
+  0:* :: ""
+} 0
+
+INT parray1d_timederivs[100] "Time derivative order for each grid function"
+{
+  0:* :: ""
+} 0
+
+REAL parray1d_x0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+REAL parray1d_y0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+REAL parray1d_z0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+INT parray1d_npoints_i "Number of grid points for the 1D grid parrays in the i direction"
+{
+  0:100 :: ""
+} 10
+
+REAL parray1d_dx_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray1d_dy_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray1d_dz_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+
+
+INT nparrays2d "Number of 2D parallel grid arrays"
+{
+  0:100 :: ""
+} 0
+
+STRING parray2d_vars[100] "Names of the grid functions that should be interpolated on a 2D grid" STEERABLE=always
+{
+  "^$"                                                 :: "do not interpolate"
+  "^[A-Za-z][A-Za-z0-9_]*[:][:][A-Za-z][A-Za-z0-9_]*$" :: "grid function name"
+} ""
+
+REAL parray2d_x0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+REAL parray2d_y0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+REAL parray2d_z0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+INT parray2d_npoints_i "Number of grid points for the 2D parallel grid arrays in the i direction"
+{
+  0:100 :: ""
+} 10
+
+REAL parray2d_dx_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray2d_dy_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray2d_dz_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+INT parray2d_npoints_j "Number of grid points for the 2D parallel grid arrays in the j direction"
+{
+  0:100 :: ""
+} 10
+
+REAL parray2d_dx_j "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray2d_dy_j "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray2d_dz_j "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+
+
+INT nparrays3d "Number of 3D grid parrays"
+{
+  0:100 :: ""
+} 0
+
+STRING parray3d_vars[100] "Names of the grid functions that should be interpolated on a 3D grid" STEERABLE=always
+{
+  "^$"                                                 :: "do not interpolate"
+  "^[A-Za-z][A-Za-z0-9_]*[:][:][A-Za-z][A-Za-z0-9_]*$" :: "grid function name"
+} ""
+
+REAL parray3d_x0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+REAL parray3d_y0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+REAL parray3d_z0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+INT parray3d_npoints_i "Number of grid points for the 3D parallel grid arrays in the i direction"
+{
+  0:100 :: ""
+} 10
+
+REAL parray3d_dx_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray3d_dy_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray3d_dz_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+INT parray3d_npoints_j "Number of grid points for the 3D parallel grid arrays in the j direction"
+{
+  0:100 :: ""
+} 10
+
+REAL parray3d_dx_j "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray3d_dy_j "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray3d_dz_j "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+INT parray3d_npoints_k "Number of grid points for the 3D parallel grid arrays in the k direction"
+{
+  0:100 :: ""
+} 10
+
+REAL parray3d_dx_k "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray3d_dy_k "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL parray3d_dz_k "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
diff --git a/schedule.ccl b/schedule.ccl
index 1ae3eb1..bafca8c 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -4,6 +4,6 @@ SCHEDULE InterpToArray AT analysis
 {
   LANG: C
   OPTIONS: global
-  STORAGE: scalars arrays1d arrays2d arrays3d
-  TRIGGERS: scalars arrays1d arrays2d arrays3d
+  STORAGE: scalars arrays1d arrays2d arrays3d parrays1d parrays2d parrays3d
+  TRIGGERS: scalars arrays1d arrays2d arrays3d parrays1d parrays2d parrays3d
 } "Interpolate to grid arrays"
diff --git a/src/interp.c b/src/interp.c
index e2dbbec..69aac42 100644
--- a/src/interp.c
+++ b/src/interp.c
@@ -29,12 +29,14 @@ InterpToArray (CCTK_ARGUMENTS)
   CCTK_INT * restrict operation_codes;
   CCTK_INT * restrict time_deriv_order;
   
+  int group;
+  cGroupDynamicData dyndata;
+  
   int nvars;
   int npoints;
   
   int n;
   int i, j, k;
-  int d;
   int ierr;
   
   
@@ -73,37 +75,38 @@ InterpToArray (CCTK_ARGUMENTS)
       coordsz[n] = scalar_z0;
       ++n;
       assert (n == npoints);
-    
+      
       inputs = malloc (nvars * sizeof * inputs);
       assert (inputs);
-    
+      
       for (n=0; n<nvars; ++n) {
         inputs[n] = CCTK_VarIndex (scalar_vars[n]);
         if (inputs[n] < 0) {
           inputs[n] = -1;
         }
       }
-    
+      
       output_types = malloc (nvars * sizeof * output_types);
       assert (output_types);
-    
+      
       for (n=0; n<nvars; ++n) {
         output_types[n] = CCTK_VARIABLE_REAL;
       }
-    
+      
       outputs = malloc (nvars * sizeof * outputs);
       assert (outputs);
-    
+      
       for (n=0; n<nvars; ++n) {
         outputs[n] = &scalars[npoints * n];
       }
-
+      
       ierr = CCTK_InterpGridArrays
         (cctkGH, 3, interpolator, options_table, coord_handle,
          npoints, CCTK_VARIABLE_REAL, coords,
          nvars, inputs,
          nvars, output_types, outputs);
-    
+      assert (! ierr);
+      
       free (coordsx);
       free (coordsy);
       free (coordsz);
@@ -140,10 +143,10 @@ InterpToArray (CCTK_ARGUMENTS)
         ++n;
       }
       assert (n == npoints);
-    
+      
       inputs = malloc (nvars * sizeof * inputs);
       assert (inputs);
-    
+      
       for (n=0; n<nvars; ++n) {
         inputs[n] = CCTK_VarIndex (array1d_vars[n]);
         assert (inputs[n] >= 0);
@@ -151,56 +154,56 @@ InterpToArray (CCTK_ARGUMENTS)
           inputs[n] = -1;
         }
       }
-    
+      
       output_types = malloc (nvars * sizeof * output_types);
       assert (output_types);
-    
+      
       for (n=0; n<nvars; ++n) {
         output_types[n] = CCTK_VARIABLE_REAL;
       }
-    
+      
       outputs = malloc (nvars * sizeof * outputs);
       assert (outputs);
-    
+      
       for (n=0; n<nvars; ++n) {
         outputs[n] = &arrays1d[npoints * n];
       }
-
+      
       operation_codes = malloc (nvars * sizeof * operation_codes);
       assert (operation_codes);
-    
+      
       for (n=0; n<nvars; ++n) {
         operation_codes[n] = array1d_spacederivs[n];
       }
-    
+      
       ierr = Util_TableSetIntArray
         (options_table, nvars, operation_codes, "operation_codes");
       assert (! ierr);
-
+      
       time_deriv_order = malloc (nvars * sizeof * time_deriv_order);
       assert (time_deriv_order);
-    
+      
       for (n=0; n<nvars; ++n) {
         time_deriv_order[n] = array1d_timederivs[n];
       }
-    
+      
       ierr = Util_TableSetIntArray
         (options_table, nvars, time_deriv_order, "time_deriv_order");
       assert (! ierr);
-
+      
       ierr = CCTK_InterpGridArrays
         (cctkGH, 3, interpolator, options_table, coord_handle,
          npoints, CCTK_VARIABLE_REAL, coords,
          nvars, inputs,
          nvars, output_types, outputs);
       assert (! ierr);
-    
+      
       ierr = Util_TableDeleteKey (options_table, "operation_codes");
       assert (! ierr);
-    
+      
       ierr = Util_TableDeleteKey (options_table, "time_deriv_order");
       assert (! ierr);
-
+      
       free (coordsx);
       free (coordsy);
       free (coordsz);
@@ -219,7 +222,7 @@ InterpToArray (CCTK_ARGUMENTS)
     nvars = narrays2d;
     if (nvars > 0) {
       npoints = array2d_npoints_i * array2d_npoints_j;
-    
+      
       coordsx = malloc (npoints * sizeof * coordsx);
       assert (coordsx);
       coordsy = malloc (npoints * sizeof * coordsy);
@@ -229,7 +232,7 @@ InterpToArray (CCTK_ARGUMENTS)
       coords[0] = coordsx;
       coords[1] = coordsy;
       coords[2] = coordsz;
-    
+      
       n = 0;
       for (j=0; j<array2d_npoints_j; ++j) {
         for (i=0; i<array2d_npoints_i; ++i) {
@@ -241,37 +244,38 @@ InterpToArray (CCTK_ARGUMENTS)
         }
       }
       assert (n == npoints);
-    
+      
       inputs = malloc (nvars * sizeof * inputs);
       assert (inputs);
-    
+      
       for (n=0; n<nvars; ++n) {
         inputs[n] = CCTK_VarIndex (array2d_vars[n]);
         if (inputs[n] < 0) {
           inputs[n] = -1;
         }
       }
-    
+      
       output_types = malloc (nvars * sizeof * output_types);
       assert (output_types);
-    
+      
       for (n=0; n<nvars; ++n) {
         output_types[n] = CCTK_VARIABLE_REAL;
       }
-    
+      
       outputs = malloc (nvars * sizeof * outputs);
       assert (outputs);
-    
+      
       for (n=0; n<nvars; ++n) {
         outputs[n] = &arrays2d[npoints * n];
       }
-    
+      
       ierr = CCTK_InterpGridArrays
         (cctkGH, 3, interpolator, options_table, coord_handle,
          npoints, CCTK_VARIABLE_REAL, coords,
          nvars, inputs,
          nvars, output_types, outputs);
-    
+      assert (! ierr);
+      
       free (coordsx);
       free (coordsy);
       free (coordsz);
@@ -288,7 +292,7 @@ InterpToArray (CCTK_ARGUMENTS)
     nvars = narrays3d;
     if (nvars > 0) {
       npoints = array3d_npoints_i * array3d_npoints_j * array3d_npoints_k;
-    
+      
       coordsx = malloc (npoints * sizeof * coordsx);
       assert (coordsx);
       coordsy = malloc (npoints * sizeof * coordsy);
@@ -298,7 +302,7 @@ InterpToArray (CCTK_ARGUMENTS)
       coords[0] = coordsx;
       coords[1] = coordsy;
       coords[2] = coordsz;
-    
+      
       n = 0;
       for (k=0; k<array3d_npoints_k; ++k) {
         for (j=0; j<array3d_npoints_j; ++j) {
@@ -312,37 +316,294 @@ InterpToArray (CCTK_ARGUMENTS)
         }
       }
       assert (n == npoints);
-    
+      
       inputs = malloc (nvars * sizeof * inputs);
       assert (inputs);
-    
+      
       for (n=0; n<nvars; ++n) {
         inputs[n] = CCTK_VarIndex (array3d_vars[n]);
         if (inputs[n] < 0) {
           inputs[n] = -1;
         }
       }
-    
+      
       output_types = malloc (nvars * sizeof * output_types);
       assert (output_types);
-    
+      
       for (n=0; n<nvars; ++n) {
         output_types[n] = CCTK_VARIABLE_REAL;
       }
-    
+      
       outputs = malloc (nvars * sizeof * outputs);
       assert (outputs);
-    
+      
       for (n=0; n<nvars; ++n) {
         outputs[n] = &arrays3d[npoints * n];
       }
-    
+      
       ierr = CCTK_InterpGridArrays
         (cctkGH, 3, interpolator, options_table, coord_handle,
          npoints, CCTK_VARIABLE_REAL, coords,
          nvars, inputs,
          nvars, output_types, outputs);
-    
+      assert (! ierr);
+      
+      free (coordsx);
+      free (coordsy);
+      free (coordsz);
+      free (inputs);
+      free (output_types);
+      free (outputs);
+    }
+  }
+  
+  
+  
+  /* Parallel 1D Arrays */
+  {
+    nvars = nparrays1d;
+    if (nvars > 0) {
+      group = CCTK_GroupIndex ("InterpToArray::parray1d");
+      assert (group >= 0);
+      ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
+      assert (! ierr);
+      assert (dyndata.dim == 1);
+      npoints = dyndata.lsh[0];
+      
+      coordsx = malloc (npoints * sizeof * coordsx);
+      assert (coordsx);
+      coordsy = malloc (npoints * sizeof * coordsy);
+      assert (coordsy);
+      coordsz = malloc (npoints * sizeof * coordsz);
+      assert (coordsz);
+      coords[0] = coordsx;
+      coords[1] = coordsy;
+      coords[2] = coordsz;
+      
+      n = 0;
+      for (i=0; i<parray1d_npoints_i; ++i) {
+        assert (n <= npoints);
+        coordsx[n] = parray1d_x0 + (dyndata.lbnd[0] + i) * parray1d_dx_i;
+        coordsy[n] = parray1d_y0 + (dyndata.lbnd[0] + i) * parray1d_dy_i;
+        coordsz[n] = parray1d_z0 + (dyndata.lbnd[0] + i) * parray1d_dz_i;
+        ++n;
+      }
+      assert (n == npoints);
+      
+      inputs = malloc (nvars * sizeof * inputs);
+      assert (inputs);
+      
+      for (n=0; n<nvars; ++n) {
+        inputs[n] = CCTK_VarIndex (parray1d_vars[n]);
+        assert (inputs[n] >= 0);
+        if (inputs[n] < 0) {
+          inputs[n] = -1;
+        }
+      }
+      
+      output_types = malloc (nvars * sizeof * output_types);
+      assert (output_types);
+      
+      for (n=0; n<nvars; ++n) {
+        output_types[n] = CCTK_VARIABLE_REAL;
+      }
+      
+      outputs = malloc (nvars * sizeof * outputs);
+      assert (outputs);
+      
+      for (n=0; n<nvars; ++n) {
+        outputs[n] = &parrays1d[npoints * n];
+      }
+      
+      operation_codes = malloc (nvars * sizeof * operation_codes);
+      assert (operation_codes);
+      
+      for (n=0; n<nvars; ++n) {
+        operation_codes[n] = parray1d_spacederivs[n];
+      }
+      
+      ierr = Util_TableSetIntArray
+        (options_table, nvars, operation_codes, "operation_codes");
+      assert (! ierr);
+      
+      time_deriv_order = malloc (nvars * sizeof * time_deriv_order);
+      assert (time_deriv_order);
+      
+      for (n=0; n<nvars; ++n) {
+        time_deriv_order[n] = parray1d_timederivs[n];
+      }
+      
+      ierr = Util_TableSetIntArray
+        (options_table, nvars, time_deriv_order, "time_deriv_order");
+      assert (! ierr);
+      
+      ierr = CCTK_InterpGridArrays
+        (cctkGH, 3, interpolator, options_table, coord_handle,
+         npoints, CCTK_VARIABLE_REAL, coords,
+         nvars, inputs,
+         nvars, output_types, outputs);
+      assert (! ierr);
+      
+      ierr = Util_TableDeleteKey (options_table, "operation_codes");
+      assert (! ierr);
+      
+      ierr = Util_TableDeleteKey (options_table, "time_deriv_order");
+      assert (! ierr);
+      
+      free (coordsx);
+      free (coordsy);
+      free (coordsz);
+      free (inputs);
+      free (output_types);
+      free (outputs);
+      free (operation_codes);
+      free (time_deriv_order);
+    }
+  }
+  
+  
+  
+  /* Parallel 2D Arrays */
+  {
+    nvars = nparrays2d;
+    if (nvars > 0) {
+      group = CCTK_GroupIndex ("InterpToArray::parray2d");
+      assert (group >= 0);
+      ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
+      assert (! ierr);
+      assert (dyndata.dim == 2);
+      npoints = dyndata.lsh[0] * dyndata.lsh[1];
+      
+      coordsx = malloc (npoints * sizeof * coordsx);
+      assert (coordsx);
+      coordsy = malloc (npoints * sizeof * coordsy);
+      assert (coordsy);
+      coordsz = malloc (npoints * sizeof * coordsz);
+      assert (coordsz);
+      coords[0] = coordsx;
+      coords[1] = coordsy;
+      coords[2] = coordsz;
+      
+      n = 0;
+      for (j=0; j<parray2d_npoints_j; ++j) {
+        for (i=0; i<parray2d_npoints_i; ++i) {
+          assert (n <= npoints);
+          coordsx[n] = parray2d_x0 + (dyndata.lbnd[0] + i) * parray2d_dx_i + (dyndata.lbnd[1] + j) * parray2d_dx_j;
+          coordsy[n] = parray2d_y0 + (dyndata.lbnd[0] + i) * parray2d_dy_i + (dyndata.lbnd[1] + j) * parray2d_dy_j;
+          coordsz[n] = parray2d_z0 + (dyndata.lbnd[0] + i) * parray2d_dz_i + (dyndata.lbnd[1] + j) * parray2d_dz_j;
+          ++n;
+        }
+      }
+      assert (n == npoints);
+      
+      inputs = malloc (nvars * sizeof * inputs);
+      assert (inputs);
+      
+      for (n=0; n<nvars; ++n) {
+        inputs[n] = CCTK_VarIndex (parray2d_vars[n]);
+        if (inputs[n] < 0) {
+          inputs[n] = -1;
+        }
+      }
+      
+      output_types = malloc (nvars * sizeof * output_types);
+      assert (output_types);
+      
+      for (n=0; n<nvars; ++n) {
+        output_types[n] = CCTK_VARIABLE_REAL;
+      }
+      
+      outputs = malloc (nvars * sizeof * outputs);
+      assert (outputs);
+      
+      for (n=0; n<nvars; ++n) {
+        outputs[n] = &parrays2d[npoints * n];
+      }
+      
+      ierr = CCTK_InterpGridArrays
+        (cctkGH, 3, interpolator, options_table, coord_handle,
+         npoints, CCTK_VARIABLE_REAL, coords,
+         nvars, inputs,
+         nvars, output_types, outputs);
+      assert (! ierr);
+      
+      free (coordsx);
+      free (coordsy);
+      free (coordsz);
+      free (inputs);
+      free (output_types);
+      free (outputs);
+    }
+  }
+  
+  
+  
+  /* Parallel 3D Arrays */
+  {
+    nvars = nparrays3d;
+    if (nvars > 0) {
+      group = CCTK_GroupIndex ("InterpToArray::parray3d");
+      assert (group >= 0);
+      ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
+      assert (! ierr);
+      assert (dyndata.dim == 3);
+      npoints = dyndata.lsh[0] * dyndata.lsh[1] * dyndata.lsh[2];
+      
+      coordsx = malloc (npoints * sizeof * coordsx);
+      assert (coordsx);
+      coordsy = malloc (npoints * sizeof * coordsy);
+      assert (coordsy);
+      coordsz = malloc (npoints * sizeof * coordsz);
+      assert (coordsz);
+      coords[0] = coordsx;
+      coords[1] = coordsy;
+      coords[2] = coordsz;
+      
+      n = 0;
+      for (k=0; k<parray3d_npoints_k; ++k) {
+        for (j=0; j<parray3d_npoints_j; ++j) {
+          for (i=0; i<parray3d_npoints_i; ++i) {
+            assert (n <= npoints);
+            coordsx[n] = parray3d_x0 + (dyndata.lbnd[0] + i) * parray3d_dx_i + (dyndata.lbnd[1] + j) * parray3d_dx_j + (dyndata.lbnd[2] + k) * parray3d_dx_k;
+            coordsy[n] = parray3d_y0 + (dyndata.lbnd[0] + i) * parray3d_dy_i + (dyndata.lbnd[1] + j) * parray3d_dy_j + (dyndata.lbnd[2] + k) * parray3d_dy_k;
+            coordsz[n] = parray3d_z0 + (dyndata.lbnd[0] + i) * parray3d_dz_i + (dyndata.lbnd[1] + j) * parray3d_dz_j + (dyndata.lbnd[2] + k) * parray3d_dz_k;
+            ++n;
+          }
+        }
+      }
+      assert (n == npoints);
+      
+      inputs = malloc (nvars * sizeof * inputs);
+      assert (inputs);
+      
+      for (n=0; n<nvars; ++n) {
+        inputs[n] = CCTK_VarIndex (parray3d_vars[n]);
+        if (inputs[n] < 0) {
+          inputs[n] = -1;
+        }
+      }
+      
+      output_types = malloc (nvars * sizeof * output_types);
+      assert (output_types);
+      
+      for (n=0; n<nvars; ++n) {
+        output_types[n] = CCTK_VARIABLE_REAL;
+      }
+      
+      outputs = malloc (nvars * sizeof * outputs);
+      assert (outputs);
+      
+      for (n=0; n<nvars; ++n) {
+        outputs[n] = &parrays3d[npoints * n];
+      }
+      
+      ierr = CCTK_InterpGridArrays
+        (cctkGH, 3, interpolator, options_table, coord_handle,
+         npoints, CCTK_VARIABLE_REAL, coords,
+         nvars, inputs,
+         nvars, output_types, outputs);
+      assert (! ierr);
+      
       free (coordsx);
       free (coordsy);
       free (coordsz);