Replaced all the CCTK_InterpGV() calls by calls to the new interpolation API

CCTK_InterpGridArrays(). Note that you now need to also compile in and activate
a thorn which provides local interpolation operators (eg. LocalInterp from
the CactusBase arrangement).

Also evaluate the return code from CCTK_InterpGridArrays() and print a warning
in case of an error.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@324 89daf98e-ef62-4674-b946-b8ff9de2216c

1 files changed, 52 insertions, 38 deletions

diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F
index d66c214..8005843 100644
--- a/src/AHFinder_gau.F
+++ b/src/AHFinder_gau.F
@@ -2,7 +2,7 @@
   @file      AHFinder_gau.F
   @date      October 1998
   @author    Miguel Alcubierre
-  @desc 
+  @desc
              Find gaussian curvature of surface.  The gaussian
              curvature is defined as R/2, where R is the Ricci
              scalar of the induced 2-geometry of the surface.
@@ -10,7 +10,7 @@
              As an extra "goody", this routine also integrates
              the equatorial circumference of the horizon and
              two polar circumferences at phi=0 and phi=pi/2.
-  @enddesc 
+  @enddesc
 @@*/
 
 #include "cctk.h"
@@ -32,11 +32,13 @@
 
       integer i,j,k,l,m,n,p
       integer npoints
-      integer interp_handle,coord_system_handle,sum_handle
+      integer param_table_handle,interp_handle,coord_system_handle,sum_handle
       integer ierror
 
       CCTK_INT rerror,error1,error2
 
+      character(30) options_string
+
       CCTK_REAL LEGEN
       CCTK_REAL theta,phi,xp,yp,zp,rp
       CCTK_REAL cost,sint,cosp,sinp
@@ -53,8 +55,10 @@
       CCTK_REAL, dimension(2,2,2)   :: d1a,gammad,gamma
       CCTK_REAL, dimension(2,2,2,2) :: d2a,gamma2
 
-      integer    num_arrays
-      integer,   dimension(7)   :: in_array_indices
+      CCTK_POINTER, dimension(3) :: interp_coords
+      CCTK_INT, dimension(7) :: in_array_indices
+      CCTK_POINTER, dimension(7) :: out_arrays
+      CCTK_INT, dimension(7) :: out_array_type_codes
 
       CCTK_REAL, allocatable, dimension(:,:) :: rr,xa,ya,za
       CCTK_REAL, allocatable, dimension(:,:) :: txx,tyy,tzz,txy,txz,tyz
@@ -205,7 +209,7 @@
 
       call CCTK_ReductionArrayHandle(sum_handle,"sum")
 
-      if (sum_handle.lt.0) then 
+      if (sum_handle.lt.0) then
         call CCTK_WARN(1,"Cannot get handle for sum reduction ! Forgot to activate an implementation providing reduction operators ??")
       end if
 
@@ -271,7 +275,7 @@
       g12 = zero
 
       gaussian = zero
-      
+
 
 !     *********************************
 !     ***   INITIALIZE PARAMETERS   ***
@@ -630,16 +634,16 @@
 !     Reduce the errors across processors (all processors must
 !     know about this since all will participate on the interpolation
 !     below).
-      
-      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle, 
-     .   error1,rerror,CCTK_VARIABLE_INT) 
+
+      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle,
+     .   error1,rerror,CCTK_VARIABLE_INT)
 
       if (ierror.ne.0) then
          call CCTK_WARN(1,"Reduction failed!")
       end if
       error1 = rerror
 
-      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle, 
+      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle,
      .   error2,rerror,CCTK_VARIABLE_INT)
 
       if (ierror.ne.0) then
@@ -670,30 +674,32 @@
 !     should be done in the routine AHFinder_int.F).  I need to fix
 !     this sometime soon! (October 2000).
 
-!     Interpolator and coordinate system handle.
-
+!     parameter, local interpolator, and coordinate system handle.
+      param_table_handle = -1
       interp_handle = -1
       coord_system_handle = -1
 
-      if (interpolation_order .eq. 1) then
-        call CCTK_InterpHandle (interp_handle, "first-order uniform cartesian")
-      else if (interpolation_order .eq. 2) then
-        call CCTK_InterpHandle (interp_handle, "second-order uniform cartesian")
-      else if (interpolation_order .eq. 3) then
-        call CCTK_InterpHandle (interp_handle, "third-order uniform cartesian")
+      options_string = "order = " // char(ichar('0') + interpolation_order)
+      call Util_TableCreateFromString (param_table_handle, options_string)
+      if (param_table_handle .lt. 0) then
+        call CCTK_WARN(0,"Cannot create parameter table for interpolator")
       endif
 
+      call CCTK_InterpHandle (interp_handle,"Lagrange polynomial interpolation")
       if (interp_handle .lt. 0) then
-        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
+        call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
       endif
 
       call CCTK_CoordSystemHandle (coord_system_handle, "cart3d")
-
       if (coord_system_handle .lt. 0) then
-        call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
+        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
       endif
 
-!     Interpolation.
+
+!     fill in the input/output arrays for the interpolator
+      interp_coords(1) = CCTK_PointerTo(xa)
+      interp_coords(2) = CCTK_PointerTo(ya)
+      interp_coords(3) = CCTK_PointerTo(za)
 
       call CCTK_VarIndex (in_array_indices(1), "admbase::gxx")
       call CCTK_VarIndex (in_array_indices(2), "admbase::gyy")
@@ -703,21 +709,29 @@
       call CCTK_VarIndex (in_array_indices(6), "admbase::gyz")
       call CCTK_VarIndex (in_array_indices(7), "ahfinder::ahfgauss")
 
-      num_arrays = 7
-      call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
-     .                    npoints, num_arrays, num_arrays,
-     .                    xa, ya, za,
-     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                    CCTK_VARIABLE_REAL,
-     .                    in_array_indices(1), in_array_indices(2),
-     .                    in_array_indices(3), in_array_indices(4),
-     .                    in_array_indices(5), in_array_indices(6),
-     .                    in_array_indices(7),
-     .                    txx, tyy, tzz, txy, txz, tyz, gaussian,
-     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL)
+      out_arrays(1) = CCTK_PointerTo(txx)
+      out_arrays(2) = CCTK_PointerTo(tyy)
+      out_arrays(3) = CCTK_PointerTo(tzz)
+      out_arrays(4) = CCTK_PointerTo(txy)
+      out_arrays(5) = CCTK_PointerTo(txz)
+      out_arrays(6) = CCTK_PointerTo(tyz)
+      out_arrays(7) = CCTK_PointerTo(gaussian)
+
+      out_array_type_codes = CCTK_VARIABLE_REAL
+
+
+!     Interpolation.
+      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
+     .                            param_table_handle, coord_system_handle,
+     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
+     .                            7, in_array_indices,
+     .                            7, out_array_type_codes, out_arrays)
+      if (ierror < 0) then
+        call CCTK_WARN (1, "interpolator call returned an error code");
+      endif
+
+!     release parameter table
+      call Util_TableDestroy (ierror, param_table_handle)
 
 
 !     *********************************