1 files changed, 200 insertions, 0 deletions
diff --git a/Carpet/CarpetInterp2/src/fasterp.hh b/Carpet/CarpetInterp2/src/fasterp.hh
new file mode 100644
index 000000000..d2a13a35f
--- /dev/null
+++ b/Carpet/CarpetInterp2/src/fasterp.hh
@@ -0,0 +1,200 @@
+#ifndef FASTERP_HH
+#define FASTERP_HH
+
+#include <cstdlib>
+#include <vector>
+
+#include <cctk.h>
+
+#include <defs.hh>
+#include <dist.hh>
+#include <typeprops.hh>
+#include <vect.hh>
+
+
+
+namespace CarpetInterp2 {
+  
+  using namespace std;
+  
+  
+  
+  int const dim = 3;
+  
+  // An interpolation point descriptor requires (3 * (max_order+1) +
+  // 1) double precision values of memory
+  int const max_order = 5;
+  
+  
+  
+  // Map C structures to MPI datatypes
+  struct mpi_struct_descr_t {
+    int          blocklength;
+    MPI_Aint     displacement;
+    MPI_Datatype type;
+  };
+  
+  void
+  create_mpi_datatype (size_t count,
+                       mpi_struct_descr_t const descr[],
+                       MPI_Datatype & newtype);
+  
+  
+  
+  // A global location, given by its global coordinates
+  struct fasterp_glocs_t {
+    vector<CCTK_REAL> coords[dim];
+    fasterp_glocs_t (size_t const n)
+    {
+      for (int d=0; d<dim; ++d) {
+        coords[d].resize(n);
+      }
+    }
+    size_t size () const { return coords[0].size(); }
+  };
+  
+  // A local location, given by map and local coordinates
+  struct fasterp_llocs_t {
+    vector<int> maps;
+    vector<CCTK_REAL> coords[dim];
+    fasterp_llocs_t (size_t const n)
+    {
+      maps.resize(n);
+      for (int d=0; d<dim; ++d) {
+        coords[d].resize(n);
+      }
+    }
+    size_t size () const { return maps.size(); }
+  };
+  
+  // An integer location, given by map, refinement level, and
+  // component
+  struct fasterp_iloc_t {
+    int m, rl, c;
+    
+    // ivect idx;
+    int ind3d;
+    rvect offset;               // in terms of grid points
+    
+    static MPI_Datatype mpi_datatype ();
+  };
+  
+  
+  
+  struct fasterp_dest_loc_t {
+    // ivect idx;
+    int ind3d;                  // destination grid point index
+  };
+  
+  class fasterp_src_loc_t {
+    CCTK_REAL coeffs[dim][max_order+1]; // interpolation coefficients
+    bvect exact;
+    
+    // ivect idx;
+    int ind3d;                  // source grid point offset
+    
+  public:
+    void
+    calc_stencil (fasterp_iloc_t const & iloc,
+                  int order);
+    void
+    interpolate (ivect const & lsh,
+                 int order,
+                 vector<CCTK_REAL const *> const & varptrs,
+                 CCTK_REAL * restrict vals)
+      const;
+    
+  private:
+    template <int O>
+    void
+    interpolate (ivect const & lsh,
+                 vector<CCTK_REAL const *> const & varptrs,
+                 CCTK_REAL * restrict vals)
+      const;
+    template <int O0, int O1, int O2>
+    void
+    interpolate (ivect const & lsh,
+                 vector<CCTK_REAL const *> const & varptrs,
+                 CCTK_REAL * restrict vals)
+      const;
+  };
+  
+  
+  
+  // A receive descriptor, describing what is received from other
+  // processors
+  struct recv_proc_t {
+    int p;                      // sending processor
+    int offset;
+    int npoints;                // total number of received points
+  };
+
+  struct recv_descr_t {
+    vector<recv_proc_t> procs;
+    vector<int> procinds;
+    int npoints;                // total number of received points
+    
+    vector<int> index;          // gather index list
+  };
+  
+  // A send descriptor; describing what to send to other processors
+  struct send_comp_t {
+    // This structure does not exist for all components -- components
+    // which are not accessed are not described, making this a sparse
+    // data structure.  The field c contains the component number.
+    vector<fasterp_src_loc_t> locs;
+    int c;                      // source component
+  };
+  
+  struct send_rl_t {
+    vector<send_comp_t> comps;
+    vector<int> compinds;
+  };
+  
+  struct send_map_t {
+    vector<send_rl_t> rls;
+  };
+
+  struct send_proc_t {
+    // This structure does not exist for all processors -- processors
+    // with which there is no communication are not described, making
+    // this a sparse data structure.  The field p contains the
+    // processor number.
+    vector<send_map_t> maps;
+    int p;                      // receiving processor
+    int offset;
+    int npoints;                // total number of sent points
+  };
+  
+  struct send_descr_t {
+    vector<send_proc_t> procs;
+    // vector<int> procinds;
+    int npoints;                // total number of sent points
+  };
+  
+  
+  
+  class fasterp_setup_t {
+    recv_descr_t recv_descr;
+    send_descr_t send_descr;
+    int order;
+    
+  public:
+    fasterp_setup_t (cGH const * restrict cctkGH,
+                     fasterp_glocs_t const & locations,
+                     int order);
+    
+    ~ fasterp_setup_t ();
+    
+    void 
+    interpolate (cGH const * restrict cctkGH,
+                 vector<int> const & varinds,
+                 vector<CCTK_REAL *> & values)
+      const;
+  };
+  
+  
+  
+} // namespace CarpetInterp2
+
+#endif  // #define FASTERP_HH