Changes to support both the new and the old interpolation API.

A local interpolator named "uniform cartesian" is now provided
under the new API.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@164 df1f8a13-aa1d-4dd4-9681-27ded5b42416

1 files changed, 42 insertions, 35 deletions

diff --git a/src/Interpolate.c b/src/Interpolate.c
index 6e9af89..4dd6067 100644
--- a/src/Interpolate.c
+++ b/src/Interpolate.c
@@ -22,8 +22,8 @@
   @author    Jonathan Thornburg
   @hdesc     Move all local-interpolation code from LocalInterp to here
   @endhistory
-  
-  @version   $Header$
+
+  @version   $Id$
   @@*/
 
 #include <math.h>
@@ -224,12 +224,12 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_Interpolate_c)
   @endvar
   @var        dims
   @vdesc      dimensions of the input arrays
-  @vtype      int[ num_dims ]
+  @vtype      CCTK_INT[ num_dims ]
   @vio        in
   @endvar
   @var        coord
-  @vdesc      coordinates to interpolate at
-  @vtype      CCTK_REAL coord[ num_dims * num_points ]
+  @vdesc      list of coordinates to interpolate at
+  @vtype      CCTK_REAL coord[ num_dims ][ num_points ]
   @vio        in
   @endvar
   @var        origin
@@ -244,7 +244,7 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_Interpolate_c)
   @endvar
   @var        in_types
   @vdesc      CCTK variable types of input arrays
-  @vtype      int in_types[ num_arrays ]
+  @vtype      CCTK_INT in_types[ num_arrays ]
   @vio        in
   @endvar
   @var        in_arrays
@@ -254,7 +254,7 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_Interpolate_c)
   @endvar
   @var        out_types
   @vdesc      CCTK variable types of output arrays
-  @vtype      int out_types[ num_arrays ]
+  @vtype      CCTK_INT out_types[ num_arrays ]
   @vio        in
   @endvar
   @var        out_arrays
@@ -273,13 +273,13 @@ int LocalInterp_Interpolate (int order,
                             int num_points,
                             int num_dims,
                             int num_arrays,
-                            const int dims[],
-                            const CCTK_REAL coord[],
+                            const CCTK_INT dims[],
+                            const CCTK_REAL *const coord[],
                             const CCTK_REAL origin[],
                             const CCTK_REAL delta[],
-                            const int in_types[],
+                            const CCTK_INT in_types[],
                             const void *const in_arrays[],
-                            const int out_types[],
+                            const CCTK_INT out_types[],
                             void *const out_arrays[])
 {
   int retval;
@@ -287,12 +287,13 @@ int LocalInterp_Interpolate (int order,
 #if 0
   int out_of_bounds;
 #endif
-  int max_dims[MAXDIM], point[MAXDIM];
+  CCTK_INT max_dims[MAXDIM], point[MAXDIM];
   CCTK_REAL delta_inv[MAXDIM];
   CCTK_REAL below[MAXDIM];
   CCTK_REAL offset[MAXDIM];
   CCTK_REAL coeff[MAXDIM][MAXORDER + 1];
 
+
   /* verify parameters and check against our restrictions */
   retval = -1;
   if (num_dims < 1)
@@ -345,7 +346,7 @@ int LocalInterp_Interpolate (int order,
      (with the remaining elements zeroed out)
      so that we can use nested loops over MAXDIM dimensions later on */
   memset (max_dims, 0, sizeof (max_dims));
-  memcpy (max_dims, dims, (num_dims - 1) * sizeof (int));
+  memcpy (max_dims, dims, (num_dims - 1) * sizeof (*max_dims));
 
   /* zero out the coefficients and set the elements with index 'order' to one
      so that we can use nested loops over MAXDIM dimensions later on */
@@ -370,7 +371,7 @@ int LocalInterp_Interpolate (int order,
     for (i = 0; i < num_dims; i++)
     {
       /* closest grid point for stencil/molecule */
-      point[i] = floor ((coord[num_dims*n + i] - origin[i]) * delta_inv[i]
+      point[i] = floor ((coord[i][n] - origin[i]) * delta_inv[i]
                         - 0.5 * (order - 1));
 
 #if 0
@@ -402,19 +403,19 @@ int LocalInterp_Interpolate (int order,
       below[i] = origin[i] + point[i] * delta[i];
 
       /* offset from that grid point, in fractions of grid points */
-      offset[i] = (coord[num_dims*n + i] - below[i]) * delta_inv[i];
+      offset[i] = (coord[i][n] - below[i]) * delta_inv[i];
     }
 
 #ifdef LOCALINTERP_VERBOSE_DEBUG
 if (n == LocalInterp_verbose_debug_n)
-	{
-	int ii;
-	printf("out_of_bounds = %d\n", out_of_bounds);
-		for (ii = 0 ; ii < num_dims ; ++ii)
-		{
-		printf("offset[%d] = %g\n", ii, (double) offset[ii]);
-		}
-	}
+        {
+        int ii;
+        printf("out_of_bounds = %d\n", out_of_bounds);
+                for (ii = 0 ; ii < num_dims ; ++ii)
+                {
+                printf("offset[%d] = %g\n", ii, (double) offset[ii]);
+                }
+        }
 #endif /* LOCALINTERP_VERBOSE_DEBUG */
 
 #if 0
@@ -425,7 +426,7 @@ if (n == LocalInterp_verbose_debug_n)
 
 
       /* put all information into a single message string for output */
-      msg = (char *) malloc (100 + num_dims*(10 + 4*20));
+      msg = malloc (100 + num_dims*(10 + 4*20));
       sprintf (msg,
                "Interpolation stencil/molecule out of bounds at grid point [%d",
                point[0]);
@@ -516,17 +517,17 @@ if (n == LocalInterp_verbose_debug_n)
 
 #ifdef LOCALINTERP_VERBOSE_DEBUG
 if (n == LocalInterp_verbose_debug_n)
-	{
-	int ii,mm;
-		for (ii = 0 ; ii < num_dims ; ++ii)
-		{
-			for (mm = 0 ; mm <= order ; ++mm)
-			{
-			printf("coeff[%d][%d] = %g\n",
-			       ii, mm, (double) coeff[ii][mm]);
-			}
-		}
-	}
+        {
+        int ii,mm;
+                for (ii = 0 ; ii < num_dims ; ++ii)
+                {
+                        for (mm = 0 ; mm <= order ; ++mm)
+                        {
+                        printf("coeff[%d][%d] = %g\n",
+                               ii, mm, (double) coeff[ii][mm]);
+                        }
+                }
+        }
 #endif /* LOCALINTERP_VERBOSE_DEBUG */
 
     /* now loop over all arrays to interpolate at the current point */
@@ -539,6 +540,12 @@ if (n == LocalInterp_verbose_debug_n)
         continue;
       }
 
+      /* skip this array if it's a query call only */
+      if (! out_arrays[a])
+      {
+        continue;
+      }
+
       /* now do the interpolation according to the array type
          we support all kinds of CCTK_REAL* and CCTK_COMPLEX* types here */
       if (in_types[a] == CCTK_VARIABLE_REAL)