1 files changed, 78 insertions, 0 deletions
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