Support interpolating to obtain the image values. This makes it

possible to use IOJpeg together with mesh refinement.

There is a new keyword parameter "gridpoints" with possible values
"hyperslab" and "interpolate".  Hyperslabbing works only on uniform
grids, whereas interpolation also works with AMR.  Additional
parameters select where in the domain the points should be
interpolated.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusIO/IOJpeg/trunk@130 eff87b29-5268-4891-90a3-a07138403961

7 files changed, 498 insertions, 92 deletions

diff --git a/par/jpeg_amr.par b/par/jpeg_amr.par
new file mode 100644
index 0000000..7d91c1b
--- /dev/null
+++ b/par/jpeg_amr.par
@@ -0,0 +1,99 @@
+!DESC "Orbiting binary sources, example of jpeg IO with AMR"
+
+ActiveThorns = "
+        Boundary
+        CartGrid3D
+        CoordBase
+        IOUtil
+        InitBase
+        LocalInterp
+        SymBase
+        Time
+        
+        HTTPD
+        HTTPDExtra
+        Socket
+        
+        GSL
+        jpeg6b
+        
+        IOJpeg
+        
+        Carpet
+        CarpetIOBasic
+        CarpetInterp
+        CarpetLib
+        CarpetReduce
+        CarpetRegrid
+        CarpetSlab
+        LoopControl
+        
+        IDScalarWaveC
+        WaveBinarySource
+        WaveToyC
+
+        AEILocalInterp
+"
+
+Cactus::cctk_run_title = "WaveToy/httpd Example"
+Cactus::terminate      = "never"
+
+InitBase::initial_data_setup_method = "init_some_levels"
+Carpet::init_fill_timelevels        = yes
+
+Carpet::max_refinement_levels   = 2
+CarpetRegrid::refinement_levels = 2
+
+Carpet::domain_from_coordbase = yes
+Grid::type                    = "coordbase"
+
+CoordBase::spacing  = "numcells"
+CoordBase::ncells_x = 70
+CoordBase::ncells_y = 70
+CoordBase::ncells_z = 70
+CoordBase::xmin     = -1.0
+CoordBase::ymin     = -1.0
+CoordBase::zmin     = -1.0
+CoordBase::xmax     = +1.0
+CoordBase::ymax     = +1.0
+CoordBase::zmax     = +1.0
+
+IOJpeg::out_vars          = "wavetoy::phi"
+IOJpeg::mode              = "remove"
+IOJpeg::gridpoints        = "interpolate"
+
+IOJpeg::out_every         =  5
+IOJpeg::refinement_factor = 10
+#IOJpeg::colormap          = "auto-old"   # does not work with mesh refinement
+IOJpeg::colormap          = "custom"
+IOJpeg::colormap_min      =  0.0 # -1.0
+IOJpeg::colormap_max      =  1.0 # +1.0
+IOJpeg::colormap_factor   = 16
+IOJpeg::out_dir           = $parfile
+
+IOJpeg::multiply_by_radius = yes
+IOJpeg::array2d_x0         = -1.0
+IOJpeg::array2d_y0         = -1.0
+IOJpeg::array2d_z0         =  0.0
+IOJpeg::array2d_npoints_i  = 101
+IOJpeg::array2d_dx_i       =  0.02
+IOJpeg::array2d_dy_i       =  0
+IOJpeg::array2d_dz_i       =  0
+IOJpeg::array2d_npoints_j  = 101
+IOJpeg::array2d_dx_j       =  0
+IOJpeg::array2d_dy_j       =  0.02
+IOJpeg::array2d_dz_j       =  0
+
+IOBasic::outInfo_every = 1
+IOBasic::outInfo_vars  = "wavetoy::phi"
+
+Time::dtfac = 0.25
+
+IDScalarWave::initial_data = "none"
+
+WaveBinarySource::binary_omega  = 26 
+WaveBinarySource::binary_charge =  0.0001
+WaveBinarySource::binary_radius =  0.25
+WaveBinarySource::binary_size   =  0.1
+
+WaveToy::bound = "radiation"
diff --git a/par/jpeg_amr.th b/par/jpeg_amr.th
new file mode 100644
index 0000000..336cfdc
--- /dev/null
+++ b/par/jpeg_amr.th
@@ -0,0 +1,32 @@
+CactusBase/Boundary
+CactusBase/CartGrid3D
+CactusBase/CoordBase
+CactusBase/IOUtil
+CactusBase/InitBase
+CactusBase/LocalInterp
+CactusBase/SymBase
+CactusBase/Time
+
+CactusConnect/HTTPD
+CactusConnect/HTTPDExtra
+CactusConnect/Socket
+
+CactusExternal/GSL
+CactusExternal/jpeg6b
+
+CactusIO/IOJpeg
+
+Carpet/Carpet
+Carpet/CarpetIOBasic
+Carpet/CarpetInterp
+Carpet/CarpetLib
+Carpet/CarpetReduce
+Carpet/CarpetRegrid
+Carpet/CarpetSlab
+Carpet/LoopControl
+
+CactusWave/IDScalarWaveC
+CactusWave/WaveBinarySource
+CactusWave/WaveToyC
+
+AEIThorns/AEILocalInterp
diff --git a/param.ccl b/param.ccl
index 15d8e89..9fdad3f 100644
--- a/param.ccl
+++ b/param.ccl
@@ -37,6 +37,12 @@ KEYWORD mode "Output mode to use" STEERABLE = ALWAYS
   "standard" :: "Generate a file for each out_every timesteps"
 } "standard"
 
+KEYWORD gridpoints "How to access grid points" STEERABLE = RECOVER
+{
+  "hyperslab"   :: "use locations of grid points"
+  "interpolate" :: "interpolate to arbitrary points"
+} "hyperslab"
+
 
 ########################
 # Specific to jpegs
@@ -58,8 +64,9 @@ INT colormap_factor "How to scale float values to rgb color" STEERABLE = ALWAYS
 
 KEYWORD colormap "How to set the colormap" STEERABLE = ALWAYS
 {
-  "auto"   :: "Set automatically using min/max of grid variables"
-  "custom" :: "Set min/max manually"
+  "auto"     :: "Set automatically using min/max of grid variables"
+  "auto-old" :: "Set automatically using min/max of grid variables, using the old reduction interface which is still used by Carpet"
+  "custom"   :: "Set min/max manually"
 } "custom"
 
 REAL colormap_min "minimum value to be mapped to colors" STEERABLE = ALWAYS
@@ -77,9 +84,9 @@ INT refinement_factor "Refine each 2D slice by a certain factor (using interpola
 } 1
 
 
-################################
-# Choosing what planes to output
-################################
+###################################################
+# Choosing what planes to output when hyperslabbing
+###################################################
 REAL out2D_yzplane_x "x-coord for 2D planes in yz" STEERABLE = RECOVER
 {
   *:*      :: "A value between [xmin, xmax]"
@@ -113,6 +120,87 @@ INT out2D_xyplane_zi "z-index (from 0) for 2D planes in xy" STEERABLE = RECOVER
 } -1
 
 
+###################################################
+# Choosing what region to output when interpolating
+###################################################
+
+STRING interpolator_name "Name of the interpolator" STEERABLE=always
+{
+  ".*" :: "must be a registered interpolator"
+} "Lagrange polynomial interpolation"
+
+STRING interpolator_options "Options for the interpolator" STEERABLE=always
+{
+  ".*" :: "must be a valid option specification"
+} "order=2"
+
+STRING interpolator_coordinates "Coordinate system" STEERABLE=always
+{
+  ".*" :: "must be a registered coordinate system"
+} "cart3d"
+
+BOOLEAN multiply_by_radius "Multiply valus by r" STEERABLE=always
+{
+} "no"
+
+
+
+REAL array2d_x0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+REAL array2d_y0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+REAL array2d_z0 "Origin" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+INT array2d_npoints_i "Number of grid points for the 2D grid arrays in the i direction"
+{
+  0:* :: ""
+} 10
+
+REAL array2d_dx_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL array2d_dy_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL array2d_dz_i "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+INT array2d_npoints_j "Number of grid points for the 2D grid arrays in the j direction"
+{
+  0:* :: ""
+} 10
+
+REAL array2d_dx_j "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL array2d_dy_j "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+REAL array2d_dz_j "Spacing" STEERABLE=always
+{
+  0.0:* :: ""
+} 0.0
+
+
 #############################################################################
 ### import IOUtil parameters
 #############################################################################
diff --git a/src/ChooseOutput.c b/src/ChooseOutput.c
index 7366fa7..04497c5 100644
--- a/src/ChooseOutput.c
+++ b/src/ChooseOutput.c
@@ -13,6 +13,7 @@
 #include <string.h>
 
 #include "cctk.h"
+#include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 #include "CactusBase/IOUtil/src/ioutil_Utils.h"
 #include "ioJpegGH.h"
@@ -41,8 +42,6 @@ CCTK_FILEVERSION(CactusIO_IOJpeg_ChooseOutput_c)
 /********************************************************************
  ********************    External Routines   ************************
  ********************************************************************/
-void IOJpeg_ChooseOutput (const cGH *GH);
-
 
 /*@@
    @routine   IOJpeg_ChooseOutput
@@ -60,13 +59,14 @@ void IOJpeg_ChooseOutput (const cGH *GH);
    @vio       in
    @endvar
  @@*/
-void IOJpeg_ChooseOutput (const cGH *GH)
+void IOJpeg_ChooseOutput (CCTK_ARGUMENTS)
 {
   int i, maxdim;
   ioJpegGH *myGH;
   int origin_index[3];
   CCTK_REAL origin_phys[3];
-  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
 
 
   GET_SLICE (out2D_xyplane_z, out_xyplane_z, origin_index[0], origin_phys[0]);
@@ -74,7 +74,7 @@ void IOJpeg_ChooseOutput (const cGH *GH)
   GET_SLICE (out2D_yzplane_x, out_yzplane_x, origin_index[2], origin_phys[2]);
 
   maxdim = CCTK_MaxDim ();
-  myGH = (ioJpegGH *) CCTK_GHExtension (GH, "IOJpeg");
+  myGH = (ioJpegGH *) CCTK_GHExtension (cctkGH, "IOJpeg");
   myGH->sp2xyz = (int **) malloc (3 * sizeof (int *));
 
   for (i = 0; i < maxdim; i++)
@@ -83,7 +83,7 @@ void IOJpeg_ChooseOutput (const cGH *GH)
 
     if (i > 0 && i < 3)
     {
-      IOUtil_2DPlanes (GH, i + 1, origin_index, origin_phys, myGH->sp2xyz[i]);
+      IOUtil_2DPlanes (cctkGH, i + 1, origin_index, origin_phys, myGH->sp2xyz[i]);
     }
   }
 }
diff --git a/src/Output.c b/src/Output.c
index 312e4ca..ea72af8 100644
--- a/src/Output.c
+++ b/src/Output.c
@@ -377,7 +377,7 @@ static void SetOutputFlag (int vindex, const char *optstring, void *arg)
       IOUtil_ParseOutputFrequency (CCTK_THORNSTRING, "IOJpeg::out_vars",
                                    myGH->stop_on_parse_errors,
                                    vindex, optstring,
-                                   &myGH->out_every[vindex]);
+                                   &myGH->out_every[vindex], NULL);
     }
   }
 }
diff --git a/src/Startup.c b/src/Startup.c
index 27566d8..c90472b 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -27,8 +27,6 @@ CCTK_FILEVERSION(CactusIO_IOJpeg_Startup_c)
 /********************************************************************
  ********************    External Routines   ************************
  ********************************************************************/
-void IOJpeg_Startup (void);
-
 
 /********************************************************************
  ********************    Internal Routines   ************************
@@ -49,13 +47,15 @@ static void *SetupGH (tFleshConfig *config, int conv_level, cGH *GH);
    @calls     CCTK_RegisterGHExtension
               CCTK_RegisterGHExtensionSetupGH
 @@*/
-void IOJpeg_Startup (void)
+int IOJpeg_Startup (void)
 {
   if (CCTK_MaxDim () >= 2)
   {
     CCTK_RegisterGHExtensionSetupGH (CCTK_RegisterGHExtension ("IOJpeg"),
                                      SetupGH);
   }
+
+  return 0;
 }
 
 
diff --git a/src/Write.c b/src/Write.c
index 3fa5c7d..4ada69b 100644
--- a/src/Write.c
+++ b/src/Write.c
@@ -9,6 +9,9 @@
  @@*/
 
 
+#include <assert.h>
+#include <limits.h>
+#include <math.h>
 #include <stdlib.h>
 #include <string.h>
 
@@ -118,103 +121,287 @@ int IOJpeg_Write (const cGH *GH, CCTK_INT vindex, const char *alias)
     min = colormap_min;
     max = colormap_max;
   }
-  else
+  else if (CCTK_Equals (colormap, "auto"))
   {
     i = CCTK_LocalArrayReductionHandle ("minimum");
-    i = CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_min);
+    CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_min);
     i = CCTK_LocalArrayReductionHandle ("maximum");
-    i = CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_max);
+    CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_max);
+  }
+  else if (CCTK_Equals (colormap, "auto-old"))
+  {
+    i = CCTK_ReductionHandle ("minimum");
+    CCTK_Reduce (GH, 0, i, 1, CCTK_VARIABLE_REAL, &min, 1, vindex);
+    i = CCTK_ReductionHandle ("maximum");
+    CCTK_Reduce (GH, 0, i, 1, CCTK_VARIABLE_REAL, &max, 1, vindex);
+  }
+  else
+  {
+    CCTK_WARN (CCTK_WARN_ABORT, "Setting for parameter IOpjeg::colormap not recognised");
   }
 
   /* get the number of slices to output */
   /* in general: maxdir = gdata.dim * (gdata.dim - 1) / 2; */
   maxdir = gdata.dim == 2 ? 1 : 3;
 
-  /* now do the actual I/O looping over all directions */
-  for (dir = 0; dir < maxdir; dir++)
+  if (CCTK_EQUALS(gridpoints, "hyperslab"))
   {
-    /* get the directions to span the hyperslab */
-    if (dir == 0)
+  
+    /* now do the actual I/O looping over all directions */
+    for (dir = 0; dir < maxdir; dir++)
     {
-      dir_i = 0; dir_j = 1;   /* xy */
-    }
-    else if (dir == 1)
-    {
-      dir_i = 0; dir_j = 2;   /* xz */
-    }
-    else
-    {
-      dir_i = 1; dir_j = 2;   /* yz */
-    }
+      /* get the directions to span the hyperslab */
+      if (dir == 0)
+      {
+        dir_i = 0; dir_j = 1;   /* xy */
+      }
+      else if (dir == 1)
+      {
+        dir_i = 0; dir_j = 2;   /* xz */
+      }
+      else
+      {
+        dir_i = 1; dir_j = 2;   /* yz */
+      }
+  
+      /* set the origin using the slice center from IOUtil */
+      memset (origin, 0, sizeof (origin));
+      if (gdata.grouptype == CCTK_GF)
+      {
+        origin[maxdir-dir-1] = myGH->sp2xyz[gdata.dim - 1][dir];
+      }
+  
+      /* set the direction vector */
+      memset (direction, 0, sizeof (direction));
+      direction[dir_i] = direction[gdata.dim + dir_j] = 1;
+  
+      mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 2,
+                                                direction, origin, extent,
+                                                NULL, -1, NULL, hsize);
+      if (mapping < 0)
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Failed to define hyperslab mapping for variable '%s'",
+                    fullname);
+        continue;
+      }
+      {
+        // Sanity check
+        // (if this fails, we requested an insane number of grid points)
+        int max = INT_MAX;
+        int d;
+        for (d=0; d<2; ++d) {
+          assert (hsize[d] >= 0 && hsize[d] <= max);
+          if (hsize[d] > 0) max /= hsize[d];
+        }
+        assert (CCTK_VarTypeSize (gdata.vartype) <= max);
+      }
+      total_hsize = hsize[0] * hsize[1] * CCTK_VarTypeSize (gdata.vartype);
+      assert (total_hsize >= 0);
+      if (total_hsize <= 0)
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Selected hyperslab has zero size for variable '%s' "
+                    "direction %d", fullname, dir);
+        Hyperslab_FreeMapping (mapping);
+        continue;
+      }
+  
+      /* allocate hyperslab buffer */
+      hdata = myproc == 0 ? malloc (total_hsize) : NULL;
+  
+      /* get the hyperslab */
+      i = Hyperslab_Get (GH, mapping, 0, vindex, 0, htype, hdata);
+      if (i)
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Failed to extract hyperslab for variable '%s'", fullname);
+      }
+  
+      /* release the mapping structure */
+      Hyperslab_FreeMapping (mapping);
+  
+      /* and dump the data to file */
+      if (myproc == 0)
+      {
+        if (i == 0)
+        {
+          /* refine the 2D slice if requested */
+          outdata = refinement_factor > 1 ? RefineData (hsize, hdata) : hdata;
+          if (! outdata)
+          {
+            CCTK_WARN (1, "Failed to refine 2D array");
+            outdata = hdata;
+          }
+          WriteData (GH, vindex, alias, gdata.dim, dir, min, max, hsize, outdata);
+          if (outdata != hdata)
+          {
+            free (outdata);
+          }
+        }
+  
+        /* clean up */
+        free (hdata);
+      }
+    } /* end of looping through xyz directions */
 
-    /* set the origin using the slice center from IOUtil */
-    memset (origin, 0, sizeof (origin));
-    if (gdata.grouptype == CCTK_GF)
+  }
+  else if (CCTK_EQUALS(gridpoints, "interpolate"))
+  {
+    
+    int interpolator;
+    int options_table;
+    int coord_handle;
+    
+    CCTK_REAL * restrict coordsx;
+    CCTK_REAL * restrict coordsy;
+    CCTK_REAL * restrict coordsz;
+    CCTK_POINTER_TO_CONST coords[3];
+    CCTK_INT * restrict inputs;
+    CCTK_INT * restrict output_types;
+    CCTK_POINTER * restrict outputs;
+    CCTK_INT * restrict operation_codes;
+    CCTK_INT * restrict time_deriv_order;
+    
+    const int nvars = 1;
+    int npoints;
+    
+    int n;
+    int i, j, k;
+    int d;
+    int ierr;
+    
+    CCTK_REAL *outdata;
+    
+    
+    
+    interpolator = CCTK_InterpHandle (interpolator_name);
+    assert (interpolator >= 0);
+    
+    options_table = Util_TableCreateFromString (interpolator_options);
+    assert (options_table >= 0);
+    
+    coord_handle = CCTK_CoordSystemHandle (interpolator_coordinates);
+    assert (coord_handle >= 0);
+    
+    
+    
     {
-      origin[maxdir-dir-1] = myGH->sp2xyz[gdata.dim - 1][dir];
+      ierr = Util_TableSetInt (options_table, 1, "want_global_mode");
+      assert (! ierr);
     }
-
-    /* set the direction vector */
-    memset (direction, 0, sizeof (direction));
-    direction[dir_i] = direction[gdata.dim + dir_j] = 1;
-
-    mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 2,
-                                              direction, origin, extent,
-                                              NULL, -1, NULL, hsize);
-    if (mapping < 0)
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Failed to define hyperslab mapping for variable '%s'",
-                  fullname);
-      continue;
+    
+    
+    
+    /* 2D Arrays */
+    npoints = myproc==0 ? array2d_npoints_i * array2d_npoints_j : 0;
+    
+    coordsx = malloc (npoints * sizeof * coordsx);
+    assert (npoints==0 || coordsx);
+    coordsy = malloc (npoints * sizeof * coordsy);
+    assert (npoints==0 || coordsy);
+    coordsz = malloc (npoints * sizeof * coordsz);
+    assert (npoints==0 || coordsz);
+    coords[0] = coordsx;
+    coords[1] = coordsy;
+    coords[2] = coordsz;
+    
+    if (myproc==0) {
+      n = 0;
+      for (j=0; j<array2d_npoints_j; ++j) {
+        for (i=0; i<array2d_npoints_i; ++i) {
+          assert (n <= npoints);
+          coordsx[n] = array2d_x0 + i * array2d_dx_i + j * array2d_dx_j;
+          coordsy[n] = array2d_y0 + i * array2d_dy_i + j * array2d_dy_j;
+          coordsz[n] = array2d_z0 + i * array2d_dz_i + j * array2d_dz_j;
+          ++n;
+        }
+      }
+      assert (n == npoints);
     }
-    total_hsize = hsize[0] * hsize[1] * CCTK_VarTypeSize (gdata.vartype);
-    if (total_hsize <= 0)
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Selected hyperslab has zero size for variable '%s' "
-                  "direction %d", fullname, dir);
-      Hyperslab_FreeMapping (mapping);
-      continue;
+    
+    inputs = malloc (nvars * sizeof * inputs);
+    assert (inputs);
+    
+    for (n=0; n<nvars; ++n) {
+      inputs[n] = vindex;
+      if (inputs[n] < 0) {
+        inputs[n] = -1;
+      }
     }
-
-    /* allocate hyperslab buffer */
-    hdata = myproc == 0 ? malloc (total_hsize) : NULL;
-
-    /* get the hyperslab */
-    i = Hyperslab_Get (GH, mapping, 0, vindex, 0, htype, hdata);
-    if (i)
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Failed to extract hyperslab for variable '%s'", fullname);
+    
+    output_types = malloc (nvars * sizeof * output_types);
+    assert (output_types);
+    
+    for (n=0; n<nvars; ++n) {
+      output_types[n] = CCTK_VARIABLE_REAL;
     }
-
-    /* release the mapping structure */
-    Hyperslab_FreeMapping (mapping);
-
-    /* and dump the data to file */
-    if (myproc == 0)
-    {
-      if (i == 0)
-      {
-        /* refine the 2D slice if requested */
-        outdata = refinement_factor > 1 ? RefineData (hsize, hdata) : hdata;
-        if (! outdata)
-        {
-          CCTK_WARN (1, "Failed to refine 2D array");
-          outdata = hdata;
-        }
-        WriteData (GH, vindex, alias, gdata.dim, dir, min, max, hsize, outdata);
-        if (outdata != hdata)
-        {
-          free (outdata);
-        }
+    
+    outputs = malloc (nvars * sizeof * outputs);
+    assert (outputs);
+    
+    outdata = malloc (npoints * sizeof * outdata);
+    assert (npoints==0 || outdata);
+    
+    for (n=0; n<nvars; ++n) {
+      outputs[n] = outdata;
+    }
+    
+    ierr = CCTK_InterpGridArrays
+      (GH, 3, interpolator, options_table, coord_handle,
+       npoints, CCTK_VARIABLE_REAL, coords,
+       nvars, inputs,
+       nvars, output_types, outputs);
+    assert (! ierr);
+    
+    
+    
+    /* Multiply values by r */
+    if (myproc==0) {
+      for (n=0; n<npoints; ++n) {
+        const CCTK_REAL x = coordsx[n];
+        const CCTK_REAL y = coordsy[n];
+        const CCTK_REAL z = coordsz[n];
+        const CCTK_REAL r = sqrt (pow (x, 2) + pow (y, 2) + pow (z, 2));
+        outdata[n] *= r;
       }
+    }
+    
+    
 
-      /* clean up */
-      free (hdata);
+    free (coordsx);
+    free (coordsy);
+    free (coordsz);
+    free (inputs);
+    free (output_types);
+    free (outputs);
+    
+    
+    
+    ierr = Util_TableDestroy (options_table);
+    assert (! ierr);
+    
+    
+    
+    if (myproc==0) {
+      const int dir = 0;
+      CCTK_INT hsize[2];
+      hsize[0] = array2d_npoints_i;
+      hsize[1] = array2d_npoints_j;
+      
+      assert (CCTK_VarTypeI(vindex) == CCTK_VARIABLE_REAL);
+      
+      WriteData (GH, vindex, alias, gdata.dim, dir,
+                 min, max, hsize, outdata);
     }
-  } /* end of looping through xyz directions */
+    
+    free (outdata);
+    
+  }
+  else
+  {
+    CCTK_WARN (CCTK_WARN_ABORT, "unknown setting for parameter IOJpeg::gridpoints");
+  }
 
   free (fullname);