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

4 files changed, 274 insertions, 87 deletions

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);