Code cleanup before moving into production mode.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGHIO/IOPanda/trunk@29 38c3d835-c875-442e-b0fe-21c19ce1d001

1 files changed, 207 insertions, 0 deletions

diff --git a/src/Write.c b/src/Write.c
new file mode 100644
index 0000000..9900d96
--- /dev/null
+++ b/src/Write.c
@@ -0,0 +1,207 @@
+/*@@
+   @file      Write.c
+   @date      01 Oct 1999
+   @author    Jonghyun Lee
+   @desc      Do the actual writing of a 3D grid array.
+   @enddesc 
+   @history
+   @hendhistory
+ @@*/
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "ioPandaGH.h"
+
+
+void IOPanda_getDumpData (const cGH *GH, int vindex, int timelevel,
+                          void **outme, int *free_outme, CCTK_INT4 *geom,
+                          int element_size);
+
+
+void IOPanda_getDumpData (const cGH *GH, int vindex, int timelevel,
+                          void **outme, int *free_outme, CCTK_INT4 *geom,
+                          int element_size) 
+{
+  DECLARE_CCTK_PARAMETERS
+  int i, myproc, do_downsample, dim;
+  ioGH *ioUtilGH;
+  CCTK_REAL4 *single_ptr;
+  CCTK_REAL *real_ptr;
+  CCTK_BYTE *char_ptr;
+  CCTK_INT *int_ptr;
+  pGExtras *extras;
+  void *data = CCTK_VarDataPtrI (GH, timelevel, vindex);
+
+  /* to make the compiler happy */
+  single_ptr = NULL;
+  real_ptr = NULL;
+  char_ptr = NULL;
+  int_ptr = NULL;
+
+  /* get the processor ID */
+  myproc = CCTK_MyProc (GH);
+
+  /* get the handle for IOUtil GH extensions */
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];
+
+  /* get the pGExtras pointer as a shortcut */
+  extras = ((pGA ***) PUGH_pGH (GH)->variables) [vindex][timelevel]->extras;
+
+  /* get the dimension of the variable */
+  dim = CCTK_GroupDimI (CCTK_GroupIndexFromVarI (vindex));
+
+  do_downsample = 0;
+  for (i = 0; i < dim; i++)
+    do_downsample |= ioUtilGH->downsample [i] > 1;
+
+  /* All the downsampling code is hard-coded for 3D data.
+     For other-dimensional variables we print a warning
+     and return the non-downsampled data. */
+  if (do_downsample && dim != 3) {
+    char *fullname = CCTK_FullName (vindex);
+
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "No downsampling for %dD variable '%s' ! Downsampling is only "
+                "supported for 3D variables.", dim, fullname);
+    free (fullname);
+    do_downsample = 0;
+  }
+
+  /* Simple case if no downsampling */
+  if (! do_downsample) {
+
+    if (ioUtilGH->out_single) {
+      single_ptr = (CCTK_REAL4 *) malloc (extras->npoints*sizeof (CCTK_REAL4));
+
+      for (i = 0; i < extras->npoints; i++)
+        single_ptr [i] = (CCTK_REAL4) ((CCTK_REAL *) data) [i];
+
+      *outme = single_ptr;
+      *free_outme = 1;
+    } else {
+      *outme = data;
+      *free_outme = 0;
+    }
+
+    /* set the bounds, the sizes, and the global shape */
+    for (i = 0; i < dim; i++) {
+      geom [i + 0*dim] = extras->lb [myproc][i];
+      geom [i + 1*dim] = extras->lnsize [i];
+      geom [i + 2*dim] = extras->nsize [i];
+    }
+
+  } else {
+
+    /* NOTE: the following downsampling code is hard-coded for 3D data */
+
+    int start [3], end [3];
+    int j, k, l;
+
+    /* downsampled global shape */
+    for (i = 0; i < 3; i++) {
+      geom [i + 6] = extras->nsize [i] / ioUtilGH->downsample [i];
+      if (extras->nsize [i] % ioUtilGH->downsample [i])
+        geom [i + 6]++;
+    }
+
+    if (verbose)
+      CCTK_VInfo (CCTK_THORNSTRING, "Downsampled sizes (%d, %d, %d) -> "
+                  "(%d, %d, %d)",
+                  extras->nsize [0], extras->nsize [1], extras->nsize [2],
+                  (int) geom [6], (int) geom [7], (int) geom [8]);
+
+    /* Now figure out the local downsampling */
+    /* The local starts are the lb modded into the downsample */
+    for (i = 0; i < 3; i++) {
+      geom [i]  = extras->lb [myproc][i] / ioUtilGH->downsample [i];
+      start [i] = geom [i] * ioUtilGH->downsample [i];
+      if (start [i] <
+          extras->lb [myproc][i] + extras->ownership [PUGH_NO_STAGGER][0][i]) {
+        start [i] += ioUtilGH->downsample [i];
+        geom [i] ++;
+      }
+      end [i]   = ((extras->lb [myproc][i] +
+                    extras->ownership [PUGH_NO_STAGGER][1][i] - 1) /
+                   ioUtilGH->downsample [i]) * ioUtilGH->downsample [i];
+      geom [i+3] = (end [i] - start [i]) / ioUtilGH->downsample [i] + 1;
+    }
+
+    if (verbose) {
+      CCTK_VInfo (CCTK_THORNSTRING, "Downsample ranges (%d, %d, %d) -> "
+                  "(%d, %d, %d)",
+                  start [0], start [1], start [2], end [0], end [1], end [2]);
+      CCTK_VInfo (CCTK_THORNSTRING, "Local size/bound  (%d, %d, %d)  "
+                  "(%d, %d, %d)", (int) geom [3], (int) geom [4], (int) geom[5],
+                  (int) geom [0], (int) geom [1], (int) geom [2]);
+    }
+
+    /* compute local ranges */
+    for (i = 0; i < 3; i++) {
+      start [i] -= extras->lb [myproc][i];
+      end [i] -= extras->lb [myproc][i];
+    }
+
+    *outme = malloc (geom [3] * geom [4] * geom [5] * element_size);
+    *free_outme = 1;
+
+    /* I hate it to repeat the loops for each case label
+       but that way produces much more efficient code */
+    l = 0;
+    switch (CCTK_VarTypeI (vindex)) {
+      case CCTK_VARIABLE_CHAR:
+        char_ptr = (CCTK_BYTE *) *outme;
+        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample [2])
+          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample [1])
+            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample [0])
+              char_ptr [l++] = ((CCTK_BYTE *) data) [DATINDEX (extras, i, j, k)];
+        break;
+
+      case CCTK_VARIABLE_INT:
+        int_ptr = (CCTK_INT *) *outme;
+        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample [2])
+          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample [1])
+            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample [0])
+              int_ptr [l++] = ((CCTK_INT *) data) [DATINDEX (extras, i, j, k)];
+        break;
+
+      case CCTK_VARIABLE_REAL:
+        if (ioUtilGH->out_single)
+          single_ptr = (CCTK_REAL4 *) *outme;
+        else
+          real_ptr = (CCTK_REAL *) *outme;
+        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample [2])
+          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample [1])
+            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample [0])
+              if (ioUtilGH->out_single)
+                single_ptr [l++] = (CCTK_REAL4)
+                                  (((CCTK_REAL *) data) [DATINDEX (extras, i, j, k)]);
+              else
+                real_ptr [l++] = ((CCTK_REAL *) data) [DATINDEX (extras, i, j, k)];
+        break;
+
+      default:
+        CCTK_WARN (1, "Unsupported variable type in IOPanda_getDumpData");
+        return;
+    }
+  }
+
+#ifdef IOPANDA_DEBUG
+  printf ("Lower bound: %d", (int) geom [0]);
+  for (i = 1; i < dim; i++)
+    printf (" %d", (int) geom [i]);
+  printf ("\n");
+  printf ("Chunk size : %d", (int) geom [1*dim + 0]);
+  for (i = 1; i < dim; i++)
+    printf (" %d", (int) geom [1*dim + i]);
+  printf ("\n");
+  printf ("Global size: %d", (int) geom [2*dim + 0]);
+  for (i = 1; i < dim; i++)
+    printf (" %d", (int) geom [2*dim + i]);
+  printf ("\n");
+#endif
+}