1 files changed, 92 insertions, 2 deletions

diff --git a/src/Write.c b/src/Write.c
index a4e1612..8c416c9 100644
--- a/src/Write.c
+++ b/src/Write.c
@@ -28,6 +28,7 @@ CCTK_FILEVERSION(CactusIO_IOJpeg_Write_c)
 static void WriteData (const cGH *GH, int vindex, const char *alias, int dim,
                        int dir, CCTK_REAL min, CCTK_REAL max,
                        const CCTK_INT hsize[2], const void *hdata);
+static void *RefineData (CCTK_INT input_dims[2], const void *input_data);
 
 
 /*@@
@@ -77,7 +78,7 @@ int IOJpeg_Write (const cGH *GH, CCTK_INT vindex, const char *alias)
   char *fullname;
   int extent_int[3];
   CCTK_INT origin[3], extent[2], direction[6], hsize[2];
-  void *hdata;
+  void *hdata, *outdata;
   CCTK_REAL min, max;
   const CCTK_INT htype = CCTK_VARIABLE_REAL;
   DECLARE_CCTK_PARAMETERS
@@ -193,7 +194,18 @@ int IOJpeg_Write (const cGH *GH, CCTK_INT vindex, const char *alias)
     {
       if (i == 0)
       {
-        WriteData (GH, vindex, alias, gdata.dim, dir, min, max, hsize, hdata);
+        /* 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 */
@@ -211,6 +223,84 @@ int IOJpeg_Write (const cGH *GH, CCTK_INT vindex, const char *alias)
  ********************    Internal Routines   ************************
  ********************************************************************/
 /*@@
+   @routine    RefineData
+   @date       Thu 4 November 2004
+   @author     Thomas Radke
+   @desc
+               Refines a given 2D input array by a certain factor
+   @enddesc
+
+   @returntype void *
+   @returndesc
+               a pointer to the allocated refined output array (must be freed
+               by the caller), or NULL in case of an error
+   @endreturndesc
+ @@*/
+static void *RefineData (CCTK_INT input_dims[2], const void *input_data)
+{
+  const CCTK_REAL coord_origin[2] = {0.0, 0.0};
+  CCTK_REAL coord_delta[2];
+  const CCTK_INT array_type_codes[2] = {CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL};
+  CCTK_REAL *interp_coords[2], *refined_data;
+  int i, j, interp_handle, N_interp_points;
+  CCTK_INT output_dims[2];
+  DECLARE_CCTK_PARAMETERS
+
+
+  interp_handle = CCTK_InterpHandle ("uniform cartesian");
+  if (interp_handle < 0)
+  {
+    CCTK_WARN (1, "Couldn't get handle for interpolation operator 'uniform "
+               "cartesian'. Did you forget to activate a thorn providing "
+               "CCTK_InterpLocalUniform() ?");
+    return (NULL);
+  }
+
+  coord_delta[0] = coord_delta[1] = refinement_factor;
+
+  output_dims[0] = refinement_factor * (input_dims[0] - 1);
+  output_dims[1] = refinement_factor * (input_dims[1] - 1);
+  N_interp_points = output_dims[0] * output_dims[1];
+
+  interp_coords[0] = malloc (2 * N_interp_points * sizeof (CCTK_REAL));
+  interp_coords[1] = interp_coords[0];
+  for (i = 0; i < output_dims[0]; i++)
+  {
+    for (j = 0; j < output_dims[1]; j++)
+    {
+      *interp_coords[0]++ = i;
+      *interp_coords[1]++ = j;
+    }
+  }
+  interp_coords[0] -= N_interp_points;
+  interp_coords[1] -= N_interp_points;
+
+  refined_data = malloc (N_interp_points * sizeof (CCTK_REAL));
+  if (CCTK_InterpLocalUniform (2, interp_handle, -1,
+                               coord_origin, coord_delta, N_interp_points,
+                               CCTK_VARIABLE_REAL,
+                               (const void *const *) interp_coords,
+                               1, input_dims, array_type_codes,
+                               &input_data, 1, array_type_codes,
+                               (void *const *) &refined_data) != 0)
+  {
+    CCTK_WARN (1, "Failed to interpolate 2D array");
+    free (refined_data);
+    refined_data = NULL;
+  }
+  else
+  {
+    input_dims[0] = output_dims[0];
+    input_dims[1] = output_dims[1];
+  }
+
+  free (interp_coords[0]);
+
+  return (refined_data);
+}
+
+
+/*@@
    @routine    WriteData
    @date       Thu 18 April 2002
    @author     Thomas Radke