CarpetInterp2: Export a simple interpolator API as aliased function

3 files changed, 99 insertions, 4 deletions

diff --git a/Carpet/CarpetInterp2/interface.ccl b/Carpet/CarpetInterp2/interface.ccl
index fc822c9e1..1a8e1c177 100644
--- a/Carpet/CarpetInterp2/interface.ccl
+++ b/Carpet/CarpetInterp2/interface.ccl
@@ -41,7 +41,24 @@ CCTK_INT FUNCTION                                \
          CCTK_INT              IN  npoints,      \
          CCTK_POINTER_TO_CONST IN  globalcoords, \
          CCTK_INT        ARRAY OUT patch,        \
-         CCTK_POINTER          IN localcoords,   \
-         CCTK_POINTER          IN dadx,          \
-         CCTK_POINTER          IN ddadxdx)
+         CCTK_POINTER          IN  localcoords,  \
+         CCTK_POINTER          IN  dadx,         \
+         CCTK_POINTER          IN  ddadxdx)
 USES FUNCTION MultiPatch_GlobalToLocal
+
+
+
+CCTK_INT FUNCTION                                      \
+    InterpGridArrays                                   \
+        (CCTK_POINTER_TO_CONST IN cctkGH,              \
+         CCTK_INT              IN N_dims,              \
+         CCTK_INT              IN order,               \
+         CCTK_INT              IN N_interp_points,     \
+         CCTK_POINTER_TO_CONST IN interp_coords,       \
+         CCTK_INT              IN N_input_arrays,      \
+         CCTK_INT        ARRAY IN input_array_indices, \
+         CCTK_INT              IN N_output_arrays,     \
+         CCTK_POINTER          IN output_arrays)
+PROVIDES FUNCTION InterpGridArrays      \
+    WITH CarpetInterp2_InterpGridArrays \
+    LANGUAGE C
diff --git a/Carpet/CarpetInterp2/src/interp2.cc b/Carpet/CarpetInterp2/src/interp2.cc
new file mode 100644
index 000000000..f8a17df99
--- /dev/null
+++ b/Carpet/CarpetInterp2/src/interp2.cc
@@ -0,0 +1,78 @@
+#include <cctk.h>
+
+#include "fasterp.hh"
+
+
+
+namespace CarpetInterp2 {
+  
+  extern "C"
+  CCTK_INT
+  CarpetInterp2_InterpGridArrays (CCTK_POINTER_TO_CONST const cctkGH_,
+                                  CCTK_INT const N_dims,
+                                  CCTK_INT const order,
+                                  CCTK_INT const N_interp_points,
+                                  CCTK_POINTER_TO_CONST const interp_coords_,
+                                  CCTK_INT const N_input_arrays,
+                                  CCTK_INT const * const input_array_indices,
+                                  CCTK_INT const N_output_arrays,
+                                  CCTK_POINTER const output_arrays_)
+  {
+    // Check input values and convert types
+    cGH const * const cctkGH = static_cast<cGH const *> (cctkGH_);
+    assert (cctkGH);
+    
+    assert (N_dims == dim);
+    
+    assert (N_interp_points >= 0);
+    
+    CCTK_REAL const * const * const interp_coords =
+      static_cast<CCTK_REAL const * const *> (interp_coords_);
+    assert (interp_coords);
+    for (int d=0; d<dim; ++d) {
+      assert (N_interp_points==0 or interp_coords[d]);
+    }
+    
+    assert (N_input_arrays >= 0);
+    
+    assert (input_array_indices);
+    for (int n=0; n<N_input_arrays; ++n) {
+      assert (input_array_indices[n]>=0 and
+              input_array_indices[n]<CCTK_NumVars());
+    }
+    
+    assert (N_output_arrays >= 0);
+    assert (N_output_arrays == N_input_arrays);
+    
+    CCTK_REAL * const * const output_arrays =
+      static_cast<CCTK_REAL * const *> (output_arrays_);
+    assert (output_arrays);
+    for (int n=0; n<N_output_arrays; ++n) {
+      assert (N_output_arrays==0 or output_arrays[n]);
+    }
+    
+    // Set up interpolation
+    fasterp_glocs_t locations (N_interp_points);
+    for (int d=0; d<dim; ++d) {
+      for (int i=0; i<N_interp_points; ++i) {
+        locations.coords[d].AT(i) = interp_coords[d][i];
+      }
+    }
+    fasterp_setup_t const setup (cctkGH, locations, order);
+    
+    // Interpolate
+    vector<int> varinds (N_input_arrays);
+    for (int n=0; n<N_input_arrays; ++n) {
+      varinds.AT(n) = input_array_indices[n];
+    }
+    vector<CCTK_REAL *> values (N_input_arrays);
+    for (int n=0; n<N_input_arrays; ++n) {
+      values.AT(n) = output_arrays[n];
+    }
+    setup.interpolate (cctkGH, varinds, values);
+    
+    // Done
+    return 0;
+  }
+  
+} // namespace CarpetInterp2
diff --git a/Carpet/CarpetInterp2/src/make.code.defn b/Carpet/CarpetInterp2/src/make.code.defn
index d95df47d6..323b3652b 100644
--- a/Carpet/CarpetInterp2/src/make.code.defn
+++ b/Carpet/CarpetInterp2/src/make.code.defn
@@ -1,7 +1,7 @@
 # Main make.code.defn file for thorn CarpetInterp2
 
 # Source files in this directory
-SRCS = fasterp.cc
+SRCS = fasterp.cc interp2.cc
 
 # Subdirectories containing source files
 SUBDIRS =