Mercurial > hg > kwantix

changeset 17:72412107029b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
On the final stretch, implementing the interface in set_bdry_conds
author Jordi Gutiérrez Hermoso <jordigh@gmail.com>
date Tue, 26 Aug 2008 10:48:19 -0500
parents 29a7b95c2805
children 382e3b8e3f88
files include/interpolator.hpp include/linalg.hpp interpolator.cpp main.cpp
diffstat 4 files changed, 161 insertions(+), 122 deletions(-) [+]
line wrap: on
line diff
--- a/include/interpolator.hpp
+++ b/include/interpolator.hpp
@@ -24,7 +24,11 @@
   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
    *  \Xi} \lambda_\xi \varphi_\xi(x) \f$, where the \f$ \varphi_\xi
    *  \f$ are RBFs.
@@ -151,6 +155,15 @@
     void set_g(const map<point, double>& g);
     //@}
 
+    /** @name Reinterpolate again 
+     *
+     * As above group of functions, but for purely interpolating
+     * interpolators 
+     */
+    //@{
+    void set_bdry_values(const bdry_values& b_new);
+    //@}
+
     /** @name Linear arithmetic operators
      *
      *  These functions return a new interpolator. They are pointwise
@@ -167,6 +180,11 @@
     /// Returns a scaled interpolator.
     interpolator<RBF> operator/(double a) const;
     //@}
+
+    /** \brief Provide a unit normals object associated to the
+     *  interpolator.
+     */
+    normals get_normals() const;
     
   private:
     ///Once the matrix is defined, this function inverts it.
@@ -260,24 +278,30 @@
     interp_values operator*(const double a) const;
     interp_values operator/(const double a) const;
     //@}
-  private:
-    
+
+  protected:
     /// 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) {};
+    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;
-    ///Exception builder
-    badArgument different_rbfs(int line, string file) const;
+    ///Number of interior points
+    size_t n;
 
     template<typename RBF> friend class interpolator;
     
     /*! \name Arithmetic operators
+     *
      * For comfortable syntax.
      */
     //@{
@@ -287,8 +311,15 @@
     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>
@@ -296,6 +327,37 @@
   {
     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.
@@ -306,6 +368,7 @@
   interp_values operator*(double a, const interp_values& v); 
   interp_values operator/(double a, const interp_values& v); 
   //@}
+  
 }
 
 #endif
--- a/include/linalg.hpp
+++ b/include/linalg.hpp
@@ -396,7 +396,8 @@
 
 namespace linalg{ //Non-member functions.
   
-  /// I/O
+  /*! @name I/O
+   */
   //@{
   ///Stream insertion operator
   std::ostream& operator<<(std::ostream& os, const vector &v);
@@ -408,7 +409,10 @@
   matrix operator>>(std::istream& is, matrix& v);
   //@}
 
-  ///Some arithmetic functions for comfortable syntax.
+  /*! @name Arithmetic 
+   *
+   * Some arithmetic functions for comfortable syntax.
+   */
   //@{
 
   /// Scale a vector.
--- a/interpolator.cpp
+++ b/interpolator.cpp
@@ -441,7 +441,7 @@
       throw not_precomputed(__LINE__, __FILE__);
      
     std::vector<size_t> alpha;
-    return interp_values(precomp_values_vec[alpha],rbfs_hash);
+    return interp_values(precomp_values_vec[alpha],rbfs_hash,n);
   }
    
   template<typename RBF>
@@ -451,7 +451,7 @@
       throw not_precomputed(__LINE__, __FILE__);
     
     std::vector<size_t> alpha(k); alpha[k-1]++;
-    return interp_values(precomp_values_vec[alpha],rbfs_hash);
+    return interp_values(precomp_values_vec[alpha],rbfs_hash,n);
   }
 
   template<typename RBF>
@@ -461,7 +461,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);
+    return interp_values(precomp_values_vec[alpha],rbfs_hash,n);
   }
 
   //***************** linear operations ***************************
@@ -519,6 +519,16 @@
     u.coeffs = (this -> coeffs)*(1/a);
     return u;
   }
+
+  template<typename RBF>
+  normals interpolator<RBF>::get_normals() const{
+    if(!initted){
+      throw badArgument("Interpolator not initialised, so there are "
+			"no normals!",__FILE__,__LINE__);
+    }
+    return normals(rbfs_hash, thebvp -> get_domain() ->
+		   get_normals());  
+  }
   
   //******************* data manipulation **************************
 
@@ -562,6 +572,30 @@
   }
 
   template<typename RBF>
+  void interpolator<RBF>::set_bdry_values(const bdry_values& b_new){
+    if(!initted){
+      throw not_initted(__LINE__, __FILE__);
+    }
+    std::vector<size_t> alpha; //empty, corresponds to evaluation
+    if(precomp_values_vec.find(alpha) == precomp_values_vec.end())
+      at();
+    linalg::vector rhs = precomp_values_vec[alpha];
+    slice s1(1,n), s2(n+1,n+m);
+    rhs(s2) = b_new.v;
+    coeffs = precomp_rbfs[alpha].inv(rhs);
+    //FIXME: Don't repeat this code
+    {
+      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;
+      
+      //Pointwise map evaluations we will not compute again.
+      precomp_values.clear();
+    }
+  }
+
+  template<typename RBF>
   error_handling::badArgument
   interpolator<RBF>::not_initted(int line, string file) const
   {
@@ -576,7 +610,6 @@
   template<typename RBF>
   void interpolator<RBF>::computecoeffs()
   {
-    using namespace std;
     linalg::vector rhs(n+m);
 
     //Compute the RBF coefficients
@@ -598,7 +631,7 @@
     
     //Precompute the values for the derivatives and evaluations that
     //have already been computed.
-    {
+    {				
       typename map<std::vector<size_t>, matrix>::const_iterator I;
       
       for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++){
@@ -645,7 +678,7 @@
     if(rbfs_hash != w.rbfs_hash)
       throw different_rbfs(__LINE__, __FILE__);
     
-    return interp_values(v + w.v,rbfs_hash);
+    return interp_values(v + w.v,rbfs_hash,n);
   }
 
   interp_values 
@@ -654,7 +687,7 @@
     if(rbfs_hash != w.rbfs_hash)
       throw different_rbfs(__LINE__, __FILE__);
     
-    return interp_values(v - w.v,rbfs_hash);
+    return interp_values(v - w.v,rbfs_hash,n);
   }
 
   interp_values 
@@ -662,7 +695,7 @@
   {
     if(rbfs_hash != w.rbfs_hash)
       throw different_rbfs(__LINE__, __FILE__);
-    return interp_values(v * w.v,rbfs_hash);
+    return interp_values(v * w.v,rbfs_hash,n);
   }
 
   interp_values 
@@ -671,33 +704,41 @@
     if(rbfs_hash != w.rbfs_hash)
       throw different_rbfs(__LINE__, __FILE__);
     
-    return interp_values(v / w.v,rbfs_hash);
+    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);
+    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);
+    return interp_values(v-w, rbfs_hash,n);
   }
 
   interp_values
   interp_values::operator*(const double a) const
   {
-    return interp_values(v*a,rbfs_hash);
+    return interp_values(v*a,rbfs_hash,n);
   }
   
   interp_values
   interp_values::operator/(const double a) const
   {
-    return interp_values(v/a, rbfs_hash);
+    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 ******************
@@ -710,7 +751,7 @@
   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);
+    return interp_values(w-v.v, v.rbfs_hash,v.n);
   } 
 
   interp_values operator*(double a, const interp_values&  v)
@@ -721,7 +762,7 @@
   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);
+    return interp_values(w/v.v, v.rbfs_hash,v.n);
   } 
 
   //Instantiations
--- a/main.cpp
+++ b/main.cpp
@@ -82,32 +82,6 @@
   bvp::interp_values at() const;
 };
 
-//FIXME: This iteration must work on interpolators without boundary
-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;};
@@ -121,8 +95,7 @@
 		 size_t n, string which_interp);
 
 template<typename RBF>
-void set_bdry_conds(interpolator<RBF> &u, interpolator<RBF> &v, 
-	       const boost::shared_ptr<domain> Omega);
+void set_bdry_conds(interpolator<RBF> &u, interpolator<RBF> &v);
 
 
 //********************** Main ******************************************
@@ -141,13 +114,13 @@
 					"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));
+    shared_ptr<dirichlet_op> D(new dirichlet_op);
     zero_func Z;
+    //FIXME: Turns this into a purely interpolating interpolator,
+    //without a BVP.
     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));
+      u_bvp_init(new linear_BVP2(Omega, I, D, Z, Z)),
+      v_bvp_init(new linear_BVP2(Omega, I, D, Z, Z));
       
    
     double dt = 1e-2;
@@ -185,7 +158,7 @@
     f_v.set_dt(dt);
     f_h.set_dt(dt);
 
-    //main loop
+    //main loop, timestepping with RK4.
     size_t maxiter = 1000;
     for(size_t i = 1; i <= maxiter; i++){ 
       cout << "Now on iteration #" << i << endl;
@@ -202,14 +175,14 @@
       
       cout << "k1 done!" << endl;
 
-      set_bdry_conds(k1[0],k1[1],Omega);
+      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));
  
-      set_bdry_conds(k2[0],k2[1],Omega);
+      set_bdry_conds(k2[0],k2[1]);
 
       cout << "k2 done!" << endl;
 
@@ -218,7 +191,7 @@
       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));
       
-      set_bdry_conds(k3[0],k3[1],Omega);
+      set_bdry_conds(k3[0],k3[1]);
       
       cout << "k3 done!" << endl;
 
@@ -227,7 +200,7 @@
       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]);
       
-      set_bdry_conds(k4[0],k4[1],Omega);
+      set_bdry_conds(k4[0],k4[1]);
 
       cout << "k4 done!" << endl;
      
@@ -235,9 +208,8 @@
       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);
-     
-      //Enforce boundary conditions iteratively
-      set_bdry_conds(u1,v1,Omega);
+    
+      set_bdry_conds(u1,v1);
       
       u0 = u1;
       v0 = v1;
@@ -307,32 +279,6 @@
   return   dt*(u()*h.d(1) + h()*u.d(1) +   v()*h.d(2) + h()*v.d(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, 
@@ -374,38 +320,23 @@
   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 set_bdry_conds(interpolator<RBF> & u, interpolator<RBF> & v, 
-		    const boost::shared_ptr<domain> Omega)
+void set_bdry_conds(interpolator<RBF> & u, interpolator<RBF> & v)
 {
-  size_t m = Omega -> get_boundary().size();
-  linalg::vector u_new(m), v_new(m);
-  {
-    set<point>::const_iterator I;
-    for(I = Omega -> get_boundary().begin(); 
-	I !=  Omega->get_boundary().end();
-	I++)
-    {
-      
-    }
+  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;
+    exc.reason="u and v have incompatible domains. This should never"
+      " happen.";
+    exc.line=__LINE__; exc.file=__FILE__; throw exc;
   }
-      
-	
+
+  //Project the vector (u,v) onto the perp of (N(1),N(2)) = (-N(2),N(1)).
+  bvp::bdry_values u_bdry_new = u_bdry*N(2)*N(2) - v_bdry*N(1)*N(2);
+  bvp::bdry_values v_bdry_new = v_bdry*N(1)*N(1) - u_bdry*N(1)*N(2);
+  
+  u.set_bdry_values(u_bdry_new);
+  v.set_bdry_values(v_bdry_new);
 }