# HG changeset patch
# User Jordi Gutiérrez Hermoso <jordigh@gmail.com>
# Date 1215997810 18000
# Node ID 0531194a431db5614f8f3ad37bbc9d987591fc19
# Parent  6e06eb6ec448ec8b65e107d8d57859afd061b966
Preliminary whole domain evaluations

diff --git a/include/interpolator.hpp b/include/interpolator.hpp
--- a/include/interpolator.hpp
+++ b/include/interpolator.hpp
@@ -13,6 +13,8 @@
   using std::map;
   using boost::shared_ptr;
 
+  class interp_values;
+
   template<typename RBF>
   class interpolator : public realfunc{
   public:
@@ -49,20 +51,38 @@
     //@{
     void interpolate(const map<point, double>& Xi); 
     void interpolate(shared_ptr<linear_BVP2> bvp);
+    
+    /** Interpolate given values on the domain
+     * 
+     * The given values must match the values returned by other
+     * interpolators on compatible domains.
+     */
+    void interpolate(const interp_values& values); 
     //@}
     
     /** @name Evaluations and derivatives
+     * 
+     * Full domain evaluations and derivatives are 
      */
     //@{
-    ///Evaluation
+    ///Point evaluation
     double operator()(const point& p) const;
-    ///Evaluation
+    ///Point evaluation
     double at(const point& p) const;
 
-    /// First derivative
+    ///Full domain evaluation
+    interp_values operator()() const; 
+    ///Full domain evaluation
+    interp_values at() const; 
+
+    /// Pointwise first derivative
     double d(const point& p, size_t k) const;
-    /// Second derivatives
+    /// Full domain first erivative
+    interp_values d(size_t k) const; 
+    /// Pointwise second derivatives
     double d2(const point &p, size_t k1, size_t k2) const;
+    /// Full domain second derivatives
+    interp_values d2(size_t k1, size_t k2) const; 
     //@}
 
     /** @name Precomputation of RBFs
@@ -130,7 +150,13 @@
 
     ///Is the interpolator ready for use?
     bool initted;
-    void not_initted(int line, string file) const; //Exception thrower
+
+    /** @name Exception throwers
+     */
+    //@{
+    void not_initted(int line, string file) const; 
+    void not_precomputed(int line, string file) const;
+    //@}
     
     ///Coefficients of the RBFs
     linalg::vector coeffs;
@@ -150,26 +176,67 @@
      * evaluations of the interpolator when the interpolator's domain
      * doesn't change.
      *
-     * The keys in this vector are vectors that represent the
-     * multindex for the derivative, where missing trailing entries in
-     * the vector represents zeros. Thus, an empty vector represents
+     * The keys in this map are vectors that represent the multindex
+     * for the derivative, where missing trailing entries in the
+     * vector represents zeros. Thus, an empty vector represents
      * evaluation instead of derivatives.
      */
     mutable map<std::vector<size_t>, matrix> precomp_rbfs;
     
-    /// Have the RBFs or their derivatives been precomputed?
+    /** @name  Evaluation flags
+     *Have the RBFs or their derivatives been precomputed?
+     */
+    //@{
     bool ev_precomp;
     bool d1_precomp;
     bool d2_precomp;
+    //@}
 
     /// Precomputed values using precomp_rbfs
     mutable map<std::vector<size_t>, map<point, double> >
     precomp_values;
+    /// Same, but storing the vectors.
+    mutable map<std::vector<size_t>, linalg::vector >
+    precomp_values_vec;
     /// Invert phis\coeffs and return result vector in a map.
     map<point, double> precompute_values(const matrix& phis) const;
   };
 
-  //For comfortable syntax
+  /// Interpolated values manipulated by the interpolator class.
+  class interp_values{
+  public:
+    /** @name Arithmetic with interpolated values
+     *
+     * These operators allow pointwise arithmetic with interpolated
+     * values returned by interpolators. They must work with
+     * interpolated values from compatible interpolators
+     * (i.e. interpolators with the same RBFs).
+     */
+    //@{
+    interp_values operator+(const interp_values& w) const;
+    interp_values operator-(const interp_values& w) const;
+    interp_values operator*(const interp_values& w) const;
+    interp_values operator/(const interp_values& w) const;
+    //@}
+  private:
+    
+    /// No construction!
+    interp_values();
+    /// Interpolator creates interp_values with this
+    interp_values(const linalg::vector& v_in, size_t rbfs_hash_in) 
+      : v(v_in), rbfs_hash(rbfs_hash_in) {};
+
+    ///Interpolated values go here
+    linalg::vector v;
+    ///Two interp_values can be added iff these two values matches.
+    size_t rbfs_hash;
+    ///Exception thrower
+    void different_rbfs(int line, string file) const;
+
+    template<typename RBF> friend class interpolator;
+  };
+
+  ///For comfortable syntax
   template <typename RBF>
   interpolator<RBF> operator*(double a, const interpolator<RBF>& u)
   {
diff --git a/include/linalg.hpp b/include/linalg.hpp
--- a/include/linalg.hpp
+++ b/include/linalg.hpp
@@ -266,8 +266,10 @@
     vector operator+(const vector& w) const;
     ///Vector subtration.
     vector operator-(const vector& w) const;
-    ///Dot product.
-    double operator*(const vector& w) const; 
+    ///Element-by-element multiplication.
+    vector operator*(const vector& w) const; 
+    ///Element-by-element division.
+    vector operator/(const vector& w) const; 
     ///Computes vM where v is treated as a row vector.
     vector operator*(const matrix& M) const; 
 
@@ -284,8 +286,11 @@
     bool operator<(const vector& w) const; 
 
     //More complex operations
-   ///Euclidean norm.
+    ///Euclidean norm.
     double norm() const; 
+    
+    ///Dot product
+    double operator%(const vector& w) const;
 
     friend class vector_view;
     friend class matrix;
@@ -394,6 +399,7 @@
   vector operator*(double a, const vector& v);
   /// Euclidean norm of a vector.
   double norm(const vector& v);
+
   /// Scale a matrix.
   matrix operator*(double a, const matrix& M);
   /// Matrix inverse, computed with LU factorisation.
diff --git a/interpolator.cpp b/interpolator.cpp
--- a/interpolator.cpp
+++ b/interpolator.cpp
@@ -133,6 +133,36 @@
   }
 
   template<typename RBF>
+  void interpolator<RBF>::interpolate(const interp_values& values)
+  {
+    if(!initted){
+      badArgument exc;
+      exc.reason = "Cannot interpolate with other interpolator's "
+	"values without defining the domain.";
+      exc.line = __LINE__;
+      exc.file = __FILE__;
+      throw exc;
+    }
+    if(rbfs_hash != values.rbfs_hash){
+      badArgument exc;
+      exc.reason = "Cannot interpolate with values from incompatible "
+	"interpolator." ;
+      exc.line = __LINE__;
+      exc.file = __FILE__;
+      throw exc;
+    }
+    coeffs = M.inv(values.v);
+   
+    //Precompute the relevant values
+    {
+      typename 
+	map<std::vector<size_t>, matrix>::const_iterator I;
+      for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++)
+	precomp_values_vec[I -> first] = (I -> second)*coeffs;
+    }
+  }
+
+  template<typename RBF>
   void interpolator<RBF>::init(shared_ptr<linear_BVP2> bvp)
   {
     thebvp = bvp;
@@ -239,7 +269,8 @@
   }
   
   template<typename RBF>
-  void interpolator<RBF>::precompute_at(){
+  void interpolator<RBF>::precompute_at()
+  {
     if(!initted){
       not_initted(__LINE__, __FILE__);
     }
@@ -269,7 +300,8 @@
   }
 
   template<typename RBF>
-  double interpolator<RBF>::at(const point& p) const{
+  double interpolator<RBF>::at(const point& p) const
+  {
     if(!initted){
       not_initted(__LINE__, __FILE__);
     }
@@ -288,8 +320,37 @@
         
     return result;
   }
+
+  template<typename RBF>
+  void interpolator<RBF>::not_precomputed(int line, string file) const
+  {
+    badArgument exc;
+    exc.reason = "Cannot do whole domain evaluations without "
+      "precomputed RBFs.";
+    exc.line = line;
+    exc.file = file;
+    throw exc;
+  }
+
+  template<typename RBF>
+  interp_values interpolator<RBF>::operator()() const
+  {
+    return at();
+  }
   
   template<typename RBF>
+  interp_values interpolator<RBF>::at() const
+  {
+    if(!ev_precomp)
+      not_precomputed(__LINE__, __FILE__);
+     
+    std::vector<size_t> alpha;
+    return interp_values(precomp_values_vec[alpha],rbfs_hash);
+  }
+  
+  
+
+  template<typename RBF>
   void interpolator<RBF>::precompute_d()
   {
     if(!initted){
@@ -341,9 +402,18 @@
     double result = 0;
     for(size_t i = 1; i <= coeffs.size(); i++)
       result += coeffs(i)*rbfs[i-1].d(p,k);
-    
+   
+    return result;
+  }
 
-    return result;
+  template<typename RBF>
+  interp_values interpolator<RBF>::d(size_t k) const
+  {
+    if(!d1_precomp)
+      not_precomputed(__LINE__, __FILE__);
+    
+    std::vector<size_t> alpha(k); alpha[k-1]++;
+    return interp_values(precomp_values_vec[alpha],rbfs_hash);
   }
 
   template<typename RBF>
@@ -399,6 +469,16 @@
     return result;
   }
 
+  template<typename RBF>
+  interp_values interpolator<RBF>::d2(size_t k1, size_t k2) const
+  {
+    if(!d2_precomp)
+      not_precomputed(__LINE__, __FILE__);
+    
+    std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++;
+    return interp_values(precomp_values_vec[alpha],rbfs_hash);
+  }
+
   //********** linear operations ***************
 
   template<typename RBF>
@@ -467,7 +547,8 @@
   }
 
   template<typename RBF>
-  void interpolator<RBF>::set_g(const realfunc& g){
+  void interpolator<RBF>::set_g(const realfunc& g)
+  {
     if(!initted){
       not_initted(__LINE__, __FILE__);
     }
@@ -476,7 +557,8 @@
   }
 
   template<typename RBF>
-  void interpolator<RBF>::set_f(const map<point, double>& f){
+  void interpolator<RBF>::set_f(const map<point, double>& f)
+  {
     if(!initted){
       not_initted(__LINE__, __FILE__);
     }
@@ -485,7 +567,8 @@
   }
 
   template<typename RBF>
-  void interpolator<RBF>::set_g(const map<point, double>& g){
+  void interpolator<RBF>::set_g(const map<point, double>& g)
+  {
     if(!initted){
       not_initted(__LINE__, __FILE__);
     }
@@ -494,7 +577,8 @@
   }
 
   template<typename RBF>
-  void interpolator<RBF>::not_initted(int line, string file) const{
+  void interpolator<RBF>::not_initted(int line, string file) const
+  {
     badArgument exc;
     exc.reason = 
       "Interpolator can't interpolate without initialisation data.";
@@ -504,7 +588,8 @@
   }
 
   template<typename RBF>
-  void interpolator<RBF>::computecoeffs(){
+  void interpolator<RBF>::computecoeffs()
+  {
     using namespace std;
     linalg::vector rhs(n+m);
 
@@ -534,6 +619,52 @@
 	precomp_values[I -> first] = precompute_values(I->second);
     }
   }
+  
+  // ************************ interp_values stuff *****************
+  void interp_values::different_rbfs(int line, string file) const
+  {
+    badArgument exc;
+    exc.reason = "Can't operate on interpolated values from "
+      "incompatible interpolators";
+    exc.line = line;
+    exc.file = file;
+    throw exc;
+  }
+
+  interp_values 
+  interp_values::operator+(const interp_values& w) const
+  {
+    if(rbfs_hash != w.rbfs_hash)
+      different_rbfs(__LINE__, __FILE__);
+    
+    return interp_values(v + w.v,rbfs_hash);
+  }
+
+  interp_values 
+  interp_values::operator-(const interp_values& w) const
+  {
+    if(rbfs_hash != w.rbfs_hash)
+      different_rbfs(__LINE__, __FILE__);
+    
+    return interp_values(v - w.v,rbfs_hash);
+  }
+
+  interp_values 
+  interp_values::operator*(const interp_values& w) const
+  {
+    if(rbfs_hash != w.rbfs_hash)
+      different_rbfs(__LINE__, __FILE__);
+    return interp_values(v * w.v,rbfs_hash);
+  }
+
+  interp_values 
+  interp_values::operator/(const interp_values& w) const
+  {
+    if(rbfs_hash != w.rbfs_hash)
+      different_rbfs(__LINE__, __FILE__);
+    
+    return interp_values(v / w.v,rbfs_hash);
+  }
 
   //Instantiations
   using namespace rbf;
diff --git a/linalg.cpp b/linalg.cpp
--- a/linalg.cpp
+++ b/linalg.cpp
@@ -523,7 +523,52 @@
     return u;    
   }
 
-  double vector::operator*(const vector& w) const{
+  vector vector::operator*(const vector& w) const
+  {
+    vector u = *this;
+    try{
+      gsl_vector_mul(u.x, w.x);
+    }
+    catch(inconformantSizes& exc){
+      exc.reason = "Can't multiply vectors of different sizes.";
+      exc.file = __FILE__;
+      exc.line = __LINE__;
+      exc.n_A = x->size;
+      exc.n_B = w.x -> size;
+      throw exc;
+    }
+    return u;
+  }
+
+  vector vector::operator/(const vector& w) const
+  {
+    vector u = *this;
+    try{
+      gsl_vector_div(u.x, w.x);
+    }
+    catch(inconformantSizes& exc){
+      exc.reason = "Can't divide vectors of different sizes.";
+      exc.file = __FILE__;
+      exc.line = __LINE__;
+      exc.n_A = x->size;
+      exc.n_B = w.x -> size;
+      throw exc;
+    }
+    return u;
+  }
+
+  vector vector::operator*(const matrix& M) const
+  {
+    return T(M)* (*this);
+  }
+  
+  double vector::norm() const
+  {
+    return gsl_blas_dnrm2(x);
+  }
+
+  double vector::operator%(const vector& w) const
+  {
     double a;
     try{
       gsl_blas_ddot(x, w.x, &a);
@@ -539,14 +584,6 @@
     return a;
   }
 
-  vector vector::operator*(const matrix& M) const{
-    return M* (*this);
-  }
-  
-  double vector::norm() const{
-    return gsl_blas_dnrm2(x);
-  }
-
   //Comparison
   bool vector::operator==(const vector& w) const{
     if(x -> size != w.x -> size){
@@ -895,31 +932,39 @@
 
   // *********** non-member arithmetic functions *******************
 
-  vector operator*(double a, const vector& v){
+  vector operator*(double a, const vector& v)
+  {
     return v*a;
   }
-  double norm(const vector& v){
+  double norm(const vector& v)
+  {
     return v.norm();
   }
-  matrix operator*(double a, const matrix& M){
+
+  matrix operator*(double a, const matrix& M)
+  {
     return M*a;
   }
-  matrix inv(const matrix& A){
+  matrix inv(const matrix& A)
+  {
     return A.inv();
   }
-  matrix T(const matrix& A){
+  matrix T(const matrix& A)
+  {
     return A.T();
   }
-  double tr(const matrix& A){
+  double tr(const matrix& A)
+  {
     return A.tr();
   }
-  double det(matrix& A){
+  double det(matrix& A)
+  {
     return A.det();
   }
 
-  double cond(matrix& A){
+  double cond(matrix& A)
+  {
     return A.cond();
   }
 
-
 }
diff --git a/main-sw.cpp b/main-sw.cpp
deleted file mode 100644
--- a/main-sw.cpp
+++ /dev/null
@@ -1,381 +0,0 @@
-#include <iostream>
-#include <fstream>
-#include <map>
-#include <set>
-#include <sstream>
-
-#include <boost/shared_ptr.hpp>
-
-#include "include/linalg.hpp"
-#include "include/rbf.hpp"
-#include "include/bvp.hpp"
-#include "include/error.hpp"
-#include "include/utils.hpp"
-#include "include/diff_op.hpp"
-#include "include/interpolator.hpp"
-#include "include/func.hpp"
-
-using namespace linalg;
-using namespace bvp;
-
-#define RBF_TYPE rbf::conical
-
-//Will solve the system
-//
-//    u_t = -u*u_x - v*u_y - g*h_x         =: g1(u,v,h)
-// 
-//    v_t = -u*v_x - v*v_y - g*h_y         =: g2(u,v,h)
-//
-//    h_t = -u*h_x - h*u_x - v*h_y - h*v_y =: g3(u,v,h)
-//
-//where g is the acceleration due to gravity, u(x,y,0) = v(x,y,0) = 0,
-//and h starts out as an undisturbed surface high with a small
-//disturbance. We'll do this by directly computing the next timestep
-//by
-//
-// u_{t+1} = u_t + /\t*g1(u_t, v_t, h_t) =: f_u(u_t, v_t, h_t, /\t)
-//
-// v_{t+1} = v_t + /\t*g2(u_t, v_t, h_t) =: f_v(u_t, v_t, h_t, /\t)
-//
-// h_{t+1} = h_t + /\t*g3(u_t, v_t, h_t) =: f_h(u_t, v_t, h_t, /\t)
-//
-//Let's hope some waves propagate here. :-)
-
-const double g = 9.8; // m/s^2
-
-template<typename RBF>
-class Fgen : public realfunc{
-public:
-  void set_interps(const interpolator<RBF>& u0, 
-		   const interpolator<RBF>& v0,
-		   const interpolator<RBF>& h0);
-  void set_u(const interpolator<RBF>& u0);
-  void set_v(const interpolator<RBF>& v0);
-  void set_h(const interpolator<RBF>& h0);
-  void set_dt(double dt_in);
-protected:
-  double dt;
-  interpolator<RBF> u;
-  interpolator<RBF> v;
-  interpolator<RBF> h;
-};
-
-template<typename RBF>
-class Fu : public Fgen<RBF>{
-  using Fgen<RBF>::u;
-  using Fgen<RBF>::v;
-  using Fgen<RBF>::h;
-  using Fgen<RBF>::dt;
-public:
-  double at(const point& p) const;
-};
-
-template<typename RBF>
-class Fv : public Fgen<RBF>{
-  using Fgen<RBF>::u;
-  using Fgen<RBF>::v;
-  using Fgen<RBF>::h;
-  using Fgen<RBF>::dt;
-public:
-  double at(const point& p) const;
-};
-
-template<typename RBF>
-class Fh : public Fgen<RBF>{
-  using Fgen<RBF>::u;
-  using Fgen<RBF>::v;
-  using Fgen<RBF>::h;
-  using Fgen<RBF>::dt;
-public:
-  double at(const point& p) const;
-};
-
-class iter_neumann : public bdry_diff_op{
-public:
-  //If this is u's, set 1, for v, set 2.
-  iter_neumann(size_t u_or_v_in);
-  double at(const realfunc &f, const point &p, const vector &n) const;
-private:
-  size_t u_or_v;
-  double at(const realfunc &f, const point &p) const {return f(p);};
-};
-
-template<typename RBF>
-class Gu_or_v : public realfunc{
-public:
-  //specify the OTHER, whether u or v.
-  Gu_or_v(size_t u_or_v_in, 
-	  const shared_ptr<domain> Omega_in) : u_or_v(u_or_v_in),
-					       Omega(Omega_in) {};
-  double at(const point& p) const;
-  void set_other(const interpolator<RBF>& other_in){other = other_in;};
-private:
-  size_t u_or_v;
-  const shared_ptr<domain> Omega;
-  interpolator<RBF> other;
-};
-
-class zero_func : public realfunc{
-public:
-  double at(const point& p) const{p.size(); return 0;};
-};
-
-//************ Function declarations ****************************
-
-template<typename RBF>
-void save_interp(const shared_ptr<domain> Omega, 
-		 const interpolator<RBF>& u, 
-		 size_t n, string which_interp);
-
-template<typename RBF>
-void bdry_iter(interpolator<RBF> &u, interpolator<RBF> &v, 
-	       const boost::shared_ptr<domain> Omega);
-
-//********************** Main ******************************************
-int main(){  
-  gsl_set_error_handler(&error_handling::errorHandler);
-
-  try{
-    using namespace std;
-    using boost::shared_ptr;
-    
-    map<point, double> h_init = utils::read_pd_map("data/h_init.map");
-      
-    shared_ptr<domain> Omega(new domain("data/circ_intr.matrix",
-					"data/circ_bdry.matrix",
-					"data/circ_nrml.matrix"));
-
-    shared_ptr<Id_op> I(new Id_op);
-    shared_ptr<iter_neumann> 
-      u_bdry_op(new iter_neumann(1)), 
-      v_bdry_op(new iter_neumann(2));
-    zero_func Z;
-    shared_ptr<linear_BVP2> 
-      u_bvp_init(new linear_BVP2(Omega, I, u_bdry_op, Z, Z)),
-      v_bvp_init(new linear_BVP2(Omega, I, v_bdry_op, Z, Z));
-      
-   
-    double dt = 0.0005;
-
-    //init the interps.
-    using namespace rbf;
-    conical::set_n(5); thin_plate_spline::set_n(6); 
-    c_infty_rbf::set_epsilon(0.01);
-
-    cout << "Initialising... please wait." << endl;
-
-    interpolator<RBF_TYPE> 
-      u1, u0(u_bvp_init), 
-      v1, v0(v_bvp_init), 
-      h1, h0(h_init);
-    u1 = u0; v1 = v0; h1 = h0;
-
-    Fu<RBF_TYPE> f_u;
-    Fv<RBF_TYPE> f_v;
-    Fh<RBF_TYPE> f_h;
-
-    f_u.set_dt(dt);
-    f_v.set_dt(dt);
-    f_h.set_dt(dt);
-
-    //main loop
-    size_t maxiter = 99999;
-    for(size_t i = 1; i <= maxiter; i++){ 
-      cout << "Now on iteration #" << i << endl;
-      save_interp(Omega,u0,i,"u");
-      save_interp(Omega,v0,i,"v");
-      save_interp(Omega,h0,i,"h");
-
-      f_h.set_interps(u0,v0,h0);
-      f_u.set_interps(u0,v0,h0);
-      f_v.set_interps(u0,v0,h0);
-
-      h1.set_f(f_h);      
-      u1.set_f(f_u);
-      v1.set_f(f_v);
-
-      //Enforce boundary conditions iteratively
-      bdry_iter(u1,v1,Omega);
-      
-      u0 = u1;
-      v0 = v1;
-      h0 = h1;
-    }
-
-    return 0;
-  }
-  catch(error exc){
-    utils::show_exception(exc);
-    return 1;
-  }
-
-}
-
-//******************* Function definitions ***************************
-
-template<typename RBF>
-void Fgen<RBF>::set_interps(const interpolator<RBF>& u0, 
-			    const interpolator<RBF>& v0,
-			    const interpolator<RBF>& h0){
-  set_u(u0);
-  set_v(v0);
-  set_h(h0);
-}
-
-template<typename RBF>
-void Fgen<RBF>::set_u(const interpolator<RBF>& u0){
-  u = u0;
-}
-
-template<typename RBF>
-void Fgen<RBF>::set_v(const interpolator<RBF>& v0){
-  v = v0;
-}
-
-template<typename RBF>
-void Fgen<RBF>::set_h(const interpolator<RBF>& h0){
-  h = h0;
-}
-
-template<typename RBF>
-void Fgen<RBF>::set_dt(double dt_in){
-  dt = dt_in;
-}
-
-
-template<typename RBF>
-double Fu<RBF>::at(const point& p) const{
-  return u(p) - dt*( 
-		    u(p)*u.d(p,1) +
-		    v(p)*u.d(p,2) +
-		    g*h.d(p,1)
-		    );
-}
-
-template<typename RBF>
-double Fv<RBF>::at(const point& p) const{
-  return v(p) - dt*( 
-		    u(p)*v.d(p,1) +
-		    v(p)*v.d(p,2) +
-		    g*h.d(p,2)
-		    );
-}
-
-template<typename RBF>
-double Fh<RBF>::at(const point& p) const{
-  return h(p) - dt*(
-		    u(p)*h.d(p,1) +
-		    h(p)*u.d(p,1) +
-		 
-		    v(p)*h.d(p,2) +
-		    h(p)*v.d(p,2)
-		    );
-}
-
-iter_neumann::iter_neumann(size_t u_or_v_in)
-{
-  if(u_or_v_in < 1 or u_or_v_in > 2){
-    error_handling::badArgument exc;
-    exc.reason = "Argument to iter_neumann constructor must be 1 or 2";
-    exc.line = __LINE__;
-    exc.file = __FILE__;
-    throw exc;
-  }
-  u_or_v = u_or_v_in;
-}
-
-double iter_neumann::at(const realfunc &f, const point &p, 
-			const vector &n) const
-{
-  return n(u_or_v)*f(p);
-}
-
-template<typename RBF>
-double Gu_or_v<RBF>::at(const point& p) const
-{
-  point n = Omega -> get_normals().at(p);
-  return -n(u_or_v)*other(p);
-}
-
-
-template<typename RBF>
-void save_interp(const shared_ptr<domain> Omega, 
-		 const interpolator<RBF>& u, 
-		 size_t n, string which_interp)
-{
-  using namespace std;
-  string filename = "results/" + which_interp;
-  if(n < 10000)
-    filename += "0";
-  if(n < 1000)
-    filename += "0";
-  if(n < 100)
-    filename += "0";
-  if(n < 10)
-    filename += "0";
-  stringstream ss;
-  string n_string;
-  ss << n;
-  ss >> n_string;
-  filename += n_string + ".map";
-  ofstream ofs(filename.c_str());
-
-  slice s(1,2);
-  matrix M(Omega -> get_interior().size() +
-	   Omega -> get_boundary().size(), 3);
-  size_t k = 1;
-  for(set<point>::const_iterator i = Omega -> get_interior().begin();
-      i != Omega -> get_interior().end(); i++){
-    M(k,s) = *i;
-    M(k,3) = u(*i);
-    k++;
-  }
-  for(set<point>::const_iterator i = Omega -> get_boundary().begin();
-      i != Omega -> get_boundary().end(); i++){
-    M(k,s) = *i;
-    M(k,3) = u(*i);
-    k++;
-  }
-  ofs << M;
-}
-
-template<typename RBF>
-linalg::vector at_bdry(const interpolator<RBF>& u,
-		       const boost::shared_ptr<domain> Omega)
-{
-  using namespace linalg;
-
-  vector out(Omega -> get_boundary().size());
-  std::set<point>::const_iterator I = Omega -> get_boundary().begin();
-
-  for(size_t i = 1; i <= out.size(); i++){
-    out(i) = u(*I);
-    I++;
-  }
-
-  return out;
-}
-
-template<typename RBF>
-void bdry_iter(interpolator<RBF> &u, interpolator<RBF> &v, 
-	       const boost::shared_ptr<domain> Omega)
-{
-  Gu_or_v<RBF_TYPE> gu(2, Omega), gv(1, Omega);
-  double err = 1;
-  do{    
-    linalg::vector u_old, u_new, v_old, v_new;
-    u_old = at_bdry(u,Omega);
-    gu.set_other(v);
-    u.set_g(gu);
-    u_new = at_bdry(u,Omega);
-
-    v_old = at_bdry(v,Omega);
-    gv.set_other(u);
-    v.set_g(gv);
-    v_new = at_bdry(v,Omega);
-    
-    double err_u = norm(u_new - u_old);
-    double err_v = norm(v_new - v_old);
-    err = (err_u > err_v? err_u : err_v);
-  }while(err > 1e-2);
-}