Mercurial > hg > kwantix
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; }