changeset 18:382e3b8e3f88

Compiles, and runs, but doesn't produce desired output. Must investigate this further.
author Jordi Gutiérrez Hermoso <jordigh@gmail.com>
date Thu, 28 Aug 2008 17:27:37 -0500
parents 72412107029b
children 36a6f2e66313
files Makefile data/sw_eqs_init.m include/interpolator.hpp interpolator.cpp linalg.cpp main.cpp results/plotresults.m utils.cpp
diffstat 8 files changed, 56 insertions(+), 237 deletions(-) [+]
line wrap: on
line diff
--- a/Makefile
+++ b/Makefile
@@ -5,10 +5,11 @@
 	  -Wpointer-arith -Wcast-qual -Wcast-align -Wwrite-strings  \
 	  -fshort-enums -fno-common -Wfatal-errors 
 OBJECTFILES = interpolator.o linalg.o error.o utils.o rbf.o bvp.o \
-	      diff_op.o ddm.o func.o main.o 
+	      diff_op.o ddm.o func.o main.o interp_values.o
 HEADERFILES = include/linalg.hpp include/error.hpp include/utils.hpp \
 	      include/rbf.hpp include/bvp.hpp include/diff_op.hpp \
-	      include/interpolator.hpp include/ddm.hpp
+	      include/interpolator.hpp include/ddm.hpp \
+	      include/interp_values.hpp
 
 .SUFFIXES : .cpp .o
 
@@ -35,6 +36,7 @@
 interpolator.o : interpolator.cpp include/interpolator.hpp
 ddm.o: ddm.cpp include/ddm.hpp
 func.o : func.cpp include/func.hpp
+interp_values.o : interp_values.cpp include/interp_values.hpp
 
 main.o: main.cpp $(HEADERFILES)
 
--- a/data/sw_eqs_init.m
+++ b/data/sw_eqs_init.m
@@ -1,4 +1,4 @@
-n = 12;
+n = 12
 m = 10;
 r = linspace(0,1,n);
 
--- a/include/interpolator.hpp
+++ b/include/interpolator.hpp
@@ -18,15 +18,12 @@
 #include "linalg.hpp"
 #include "func.hpp"
 #include "diff_op.hpp"
+#include "interp_values.hpp"
 
 namespace bvp{
   using std::map;
   using boost::shared_ptr;
 
-  class interp_values;
-  class normals;
-  class bdry_values;
-
   //FIXME: Break this up into an ansatz class and an interpolator
   //class. 
   /*! \brief An interpolator \f$ \displaystyle u(x) = \sum_{\xi \in
@@ -222,6 +219,9 @@
 
     ///An RBF hasher
     size_t hash_value(const std::vector<RBF>& rbfs_in);
+
+    ///The normals at the domain's boundary
+    normals nrmls;
        
     /** \short Precomputed RBFs 
      * 
@@ -257,69 +257,6 @@
     map<point, double> valsvec2map(const linalg::vector& values) const;
   };
 
-  /// 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;
-
-    interp_values operator+(const double a) const;
-    interp_values operator-(const double a) const;
-    interp_values operator*(const double a) const;
-    interp_values operator/(const double a) const;
-    //@}
-
-  protected:
-    /// No construction!
-    interp_values();
-    /// Interpolator creates interp_values with this
-    interp_values(const linalg::vector& v_in, size_t rbfs_hash_in,
-    size_t n_in) 
-      : v(v_in), rbfs_hash(rbfs_hash_in), n(n_in) {};
-
-    ///Exception builder
-    badArgument different_rbfs(int line, string file) const;
-    
-  protected:
-    ///Interpolated values go here
-    linalg::vector v;
-    ///Two interp_values can be added iff these two values matches.
-    size_t rbfs_hash;
-    ///Number of interior points
-    size_t n;
-
-    template<typename RBF> friend class interpolator;
-    
-    /*! \name Arithmetic operators
-     *
-     * For comfortable syntax.
-     */
-    //@{
-    friend interp_values operator+(double a, const interp_values& v); 
-    friend interp_values operator-(double a, const interp_values& v); 
-    friend interp_values operator*(double a, const interp_values& v); 
-    friend interp_values operator/(double a, const interp_values& v); 
-    //@}
-  };
-  
-  /// For referring to values purely on the boundary of an interpolator.
-  class bdry_values : public interp_values{
-  public:
-    ///Initialise the boundary values with interp_values
-    bdry_values(const interp_values& v_in);
-  private:
-    bdry_values();
-  };
 
   /// For comfortable syntax
   template <typename RBF>
@@ -328,46 +265,6 @@
     return u*a;
   }
   
-  /// A unit normals class associated to an interpolator
-  class normals{
-  public:
-    /** \brief Return the kth coordinate of all normals along the
-    * boundary as interp values
-    */
-    interp_values operator()(size_t k) const;
-    
-    /*! \name comparisons 
-     *
-     * Normals are the same if they correspond to the same domain
-     */
-    //@{
-    bool operator!=(const normals& M) const;
-    bool operator==(const normals& M) const;
-    //@}
-  private:
-    template<typename RBF> friend class interpolator;
-    ///Only constructor available
-    normals(size_t rbfs_hash_in, map<point, linalg::vector>
-	    normals_in);
-    ///Exception builder
-    badArgument not_initted() const;
-    
-    ///A hash used to return inter_values objects
-    size_t rbfs_hash;
-    
-    ///Store normals in a matrix
-    linalg::matrix N;
-  };
-
-  /*! \name Arithmetic operators
-   * For comfortable syntax.
-   */
-  //@{
-  interp_values operator+(double a, const interp_values& v); 
-  interp_values operator-(double a, const interp_values& v); 
-  interp_values operator*(double a, const interp_values& v); 
-  interp_values operator/(double a, const interp_values& v); 
-  //@}
   
 }
 
--- a/interpolator.cpp
+++ b/interpolator.cpp
@@ -1,6 +1,9 @@
 #include <vector>
 #include <set>
 
+//debug
+#include <iostream>
+
 #include <utility>
 
 #include "include/interpolator.hpp"
@@ -75,7 +78,6 @@
     shared_ptr<linear_BVP2> new_bvp(new linear_BVP2(Omega, Id, D, f,g));
     
     interpolate(new_bvp);
-
   }
 
   template<typename RBF>
@@ -228,6 +230,9 @@
 
     rbfs_hash = hash_value(rbfs);
 
+    bvp::normals nrmls_tmp(rbfs_hash, normals, n);  
+    nrmls = nrmls_tmp;
+
     initted = true;
   }
 
@@ -441,7 +446,7 @@
       throw not_precomputed(__LINE__, __FILE__);
      
     std::vector<size_t> alpha;
-    return interp_values(precomp_values_vec[alpha],rbfs_hash,n);
+    return interp_values(precomp_values_vec[alpha],rbfs_hash,n,m);
   }
    
   template<typename RBF>
@@ -451,7 +456,7 @@
       throw not_precomputed(__LINE__, __FILE__);
     
     std::vector<size_t> alpha(k); alpha[k-1]++;
-    return interp_values(precomp_values_vec[alpha],rbfs_hash,n);
+    return interp_values(precomp_values_vec[alpha],rbfs_hash,n,m);
   }
 
   template<typename RBF>
@@ -461,7 +466,7 @@
       throw 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,n);
+    return interp_values(precomp_values_vec[alpha],rbfs_hash,n,m);
   }
 
   //***************** linear operations ***************************
@@ -526,8 +531,7 @@
       throw badArgument("Interpolator not initialised, so there are "
 			"no normals!",__FILE__,__LINE__);
     }
-    return normals(rbfs_hash, thebvp -> get_domain() ->
-		   get_normals());  
+    return nrmls;
   }
   
   //******************* data manipulation **************************
@@ -661,109 +665,7 @@
     return out;
   }
   
-  // ************************ interp_values stuff *****************
-  badArgument 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;
-    return exc;
-  }
 
-  interp_values 
-  interp_values::operator+(const interp_values& w) const
-  {
-    if(rbfs_hash != w.rbfs_hash)
-      throw different_rbfs(__LINE__, __FILE__);
-    
-    return interp_values(v + w.v,rbfs_hash,n);
-  }
-
-  interp_values 
-  interp_values::operator-(const interp_values& w) const
-  {
-    if(rbfs_hash != w.rbfs_hash)
-      throw different_rbfs(__LINE__, __FILE__);
-    
-    return interp_values(v - w.v,rbfs_hash,n);
-  }
-
-  interp_values 
-  interp_values::operator*(const interp_values& w) const
-  {
-    if(rbfs_hash != w.rbfs_hash)
-      throw different_rbfs(__LINE__, __FILE__);
-    return interp_values(v * w.v,rbfs_hash,n);
-  }
-
-  interp_values 
-  interp_values::operator/(const interp_values& w) const
-  {
-    if(rbfs_hash != w.rbfs_hash)
-      throw different_rbfs(__LINE__, __FILE__);
-    
-    return interp_values(v / w.v,rbfs_hash,n);
-  }
-
-  interp_values
-  interp_values::operator+(const double a) const
-  {
-    vector w(v.size(),a);
-    return interp_values(v+w,rbfs_hash,n);
-  }
-  
-  interp_values
-  interp_values::operator-(const double a) const
-  {
-    vector w(v.size(),a);
-    return interp_values(v-w, rbfs_hash,n);
-  }
-
-  interp_values
-  interp_values::operator*(const double a) const
-  {
-    return interp_values(v*a,rbfs_hash,n);
-  }
-  
-  interp_values
-  interp_values::operator/(const double a) const
-  {
-    return interp_values(v/a, rbfs_hash,n);
-  }
-
-  //*************** bdry_vals stuff ************************************
-
-  bdry_values::bdry_values(const interp_values& v_in) : interp_values(v_in){
-    //Lower part of the vector contains the boundary
-    linalg::slice s(n+1,v.size());
-    v = v(s);
-  }
-
-  //*************** Non-member global friend functions ******************
-
-  interp_values operator+(double a, const interp_values&  v)
-  { 
-    return v+a;
-  } 
-
-  interp_values operator-(double a, const interp_values&  v)
-  { 
-    linalg::vector w(v.v.size(),a);
-    return interp_values(w-v.v, v.rbfs_hash,v.n);
-  } 
-
-  interp_values operator*(double a, const interp_values&  v)
-  { 
-    return v*a;
-  } 
-
-  interp_values operator/(double a, const interp_values&  v)
-  { 
-    linalg::vector w(v.v.size(),a);
-    return interp_values(w/v.v, v.rbfs_hash,v.n);
-  } 
 
   //Instantiations
   using namespace rbf;
--- a/linalg.cpp
+++ b/linalg.cpp
@@ -735,6 +735,7 @@
       exc.reason = 
 	"Invalid arguments to slice::set. "
 	"Right endpoint must be greater than left endpoint.";
+      exc.line = __LINE__; exc.file = __FILE__;
       throw exc;
     }
     else if(k==0){
@@ -742,12 +743,14 @@
       exc.reason = 
 	"Invalid arguments to slice::set. "
 	"Cannot take zero stride.";
+      exc.line = __LINE__; exc.file = __FILE__;
       throw exc;
     }
     else if( k > b-a){
       badArgument exc;
       exc.reason = "Invalid arguments to slice::set. " 
 	"Step size cannot be greater than range.";
+      exc.line = __LINE__; exc.file = __FILE__;
       throw exc;
     }
     beg = a-1;
--- a/main.cpp
+++ b/main.cpp
@@ -3,7 +3,7 @@
 #include <map>
 #include <set>
 #include <sstream>
-
+ 
 #include <boost/shared_ptr.hpp>
 
 #include "include/linalg.hpp"
@@ -106,7 +106,7 @@
   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",
@@ -123,7 +123,7 @@
       v_bvp_init(new linear_BVP2(Omega, I, D, Z, Z));
       
    
-    double dt = 1e-2;
+    double dt = 1;
 
     //init the interps.
     using namespace rbf;
@@ -172,17 +172,25 @@
       f_u.set_interps(u0,v0,h0);
       f_v.set_interps(u0,v0,h0);
       f_h.set_interps(u0,v0,h0);
+
+      k1[0].interpolate(f_u.at());
+      k1[1].interpolate(f_v.at());
+      k1[2].interpolate(f_h.at());
+
+      //set_bdry_conds(k1[0],k1[1]);
       
       cout << "k1 done!" << endl;
-
-      set_bdry_conds(k1[0],k1[1]);
            
       //k2
       f_u.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
       f_v.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
       f_h.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
+
+      k2[0].interpolate(f_u.at());
+      k2[1].interpolate(f_v.at());
+      k2[2].interpolate(f_h.at());
  
-      set_bdry_conds(k2[0],k2[1]);
+      //set_bdry_conds(k2[0],k2[1]);
 
       cout << "k2 done!" << endl;
 
@@ -190,8 +198,12 @@
       f_u.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2));
       f_v.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2));
       f_h.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2));
+
+      k3[0].interpolate(f_u.at());
+      k3[1].interpolate(f_v.at());
+      k3[2].interpolate(f_h.at());
       
-      set_bdry_conds(k3[0],k3[1]);
+      //set_bdry_conds(k3[0],k3[1]);
       
       cout << "k3 done!" << endl;
 
@@ -199,17 +211,21 @@
       f_u.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]);
       f_v.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]);
       f_h.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]);
+
+      k4[0].interpolate(f_u.at());
+      k4[1].interpolate(f_v.at());
+      k4[2].interpolate(f_h.at());
       
-      set_bdry_conds(k4[0],k4[1]);
+      //set_bdry_conds(k4[0],k4[1]);
 
       cout << "k4 done!" << endl;
      
       //Grand finale
-      u1.set_f(u0 + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])/6);
-      v1.set_f(v0 + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])/6);
-      h1.set_f(h0 + (k1[2] + 2*k2[2] + 2*k3[2] + k4[2])/6);
+      u1.interpolate((u0 + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])/6).at());
+      v1.interpolate((v0 + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])/6).at());
+      h1.interpolate((h0 + (k1[2] + 2*k2[2] + 2*k3[2] + k4[2])/6).at());
     
-      set_bdry_conds(u1,v1);
+      //set_bdry_conds(u1,v1);
       
       u0 = u1;
       v0 = v1;
@@ -278,7 +294,7 @@
 {
   return   dt*(u()*h.d(1) + h()*u.d(1) +   v()*h.d(2) + h()*v.d(2));
 }
-
+ 
 template<typename RBF>
 void save_interp(const shared_ptr<domain> Omega, 
 		 const interpolator<RBF>& u, 
@@ -324,7 +340,7 @@
 template<typename RBF>
 void set_bdry_conds(interpolator<RBF> & u, interpolator<RBF> & v)
 {
-  bvp::bdry_values u_bdry = u.at(), v_bdry = v.at();
+  bvp::bdry_values u_bdry(u.at()), v_bdry(v.at());
   bvp::normals N=u.get_normals(), M=v.get_normals();
   if(N != M){
     failure exc;
--- a/results/plotresults.m
+++ b/results/plotresults.m
@@ -3,7 +3,7 @@
 [rr,thth] = meshgrid(r,th);
 [xx,yy] = pol2cart(thth,rr);
 
-for i = 10:10:6100
+for i = 1:1:1000
   name = "h";
   if(i < 10000);
     name = strcat(name, "0");
@@ -25,7 +25,7 @@
   surf(xx,yy,zz);
   xlim([-1,1]);
   ylim([-1,1]);
-  zlim([0,0.2]);
+##  zlim([0,0.2]);
   print(strcat(name,".png"),"-dpng");
   clear eval(name);
 endfor
--- a/utils.cpp
+++ b/utils.cpp
@@ -100,9 +100,8 @@
   void show_exception(error_handling::error exc){
     using namespace std;
 
-    cerr << "Caught an exception!" << endl;
-    cerr << exc.reason << endl;
-    cerr << "On line " << exc.line << endl;
+    cerr << exc.file << ": ""Caught an exception!" << endl;
+    cerr << exc.file << ":" << exc.line << " error: "<< exc.reason << endl;
     cerr << "From file " << exc.file << endl;
   }