# HG changeset patch # User Jordi Gutiérrez Hermoso # Date 1219963087 18000 # Node ID 36a6f2e66313a966f0d0cfb9280548c9b96543d8 # Parent 382e3b8e3f881ac36063e7791c53561f59269bcd Oops, forgot to add a new file diff --git a/include/interp_values.hpp b/include/interp_values.hpp new file mode 100644 --- /dev/null +++ b/include/interp_values.hpp @@ -0,0 +1,144 @@ +#ifndef __INTERP_VALUES_HPP__ +#define __INTERP_VALUES_HPP__ + +#include "linalg.hpp" +#include + +namespace bvp{ + using std::map; + using namespace linalg; + + class bdry_values; + + //FIXME: Give this a better name + /// 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 + */ + bdry_values operator()(size_t k) const; + + /*! \name comparisons + * + * Normals are the same if they correspond to the same domain + *(i.e. if their hash coincides) + */ + //@{ + bool operator!=(const normals& M) const; + bool operator==(const normals& M) const; + //@} + private: + template friend class interpolator; + ///Only constructor available for interpolator class + normals(size_t rbfs_hash_in, const map& + normals_in, size_t n_in); + ///Does nothing + normals(){}; + ///Exception builder + badArgument not_initted() const; + + ///A hash used to return bdry_values objects + size_t rbfs_hash; + + ///Number of interior points in the matrix, which interp_values need + size_t n; + ///Number of boundary points + size_t m; + + ///Store normals in a matrix + linalg::matrix N; + }; + + + + + /// 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; + //@} + + /*! \name Non-member 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); + //@} + + protected: + /// Interpolator creates interp_values with this + interp_values(const linalg::vector& v_in, size_t rbfs_hash_in, + size_t n_in, size_t m_in) + : v(v_in), rbfs_hash(rbfs_hash_in), n(n_in), m(m_in) {}; + + /// No construction! + interp_values(); + + protected: + ///Exception builder + badArgument different_rbfs(int line, string file) const; + + ///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; + ///Number of boundary points + size_t m; + + template friend class interpolator; + }; + + + + /// For referring to values purely on the boundary of an interpolator. + class bdry_values : public interp_values{ + public: + ///Empty initialisation + bdry_values(){}; + ///Initialise the boundary values with interp_values + bdry_values(const interp_values& v_in); + private: + friend class normals; + ///Same as for parent class + bdry_values(const linalg::vector& v_in, size_t rbfs_hash_in, + size_t n_in, size_t m_in) + : interp_values(v_in, rbfs_hash_in, n_in, m_in) {}; + }; + + + + /*! \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); + //@} +} + +#endif //__INTERP_VALUES_HPP__ diff --git a/interp_values.cpp b/interp_values.cpp new file mode 100644 --- /dev/null +++ b/interp_values.cpp @@ -0,0 +1,151 @@ +#include "include/interp_values.hpp" + +namespace bvp{ + + // ************************ 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,m); + } + + 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,m); + } + + 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,m); + } + + 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,m); + } + + interp_values + interp_values::operator+(const double a) const + { + vector w(v.size(),a); + return interp_values(v+w,rbfs_hash,n,m); + } + + interp_values + interp_values::operator-(const double a) const + { + vector w(v.size(),a); + return interp_values(v-w, rbfs_hash,n,m); + } + + interp_values + interp_values::operator*(const double a) const + { + return interp_values(v*a,rbfs_hash,n,m); + } + + interp_values + interp_values::operator/(const double a) const + { + return interp_values(v/a, rbfs_hash,n,m); + } + + //*************** bdry_vals stuff ************************************ + + bdry_values::bdry_values(const interp_values& v_in) + : interp_values(v_in) + { + if(v.size() == m+n){ + //Lower part of the vector contains the boundary + linalg::slice s(n+1,v.size()); + v = v(s); + } + else if(v.size() != m){ + throw badArgument("Cannot initialise bdry_values from given " + "interp_values", __FILE__,__LINE__); + } + } + + //*************** normals stuff ************************************* + + normals::normals(size_t rbfs_hash_in, + const map& normals_in, + size_t n_in) : rbfs_hash(rbfs_hash_in), n(n_in) + { + map::const_iterator I; + size_t rows = normals_in.size(); + size_t cols = normals_in.begin() -> second.size(); + slice s(1,cols); + size_t i = 1; + matrix N_temp(rows,cols); + for(I = normals_in.begin(); I != normals_in.end(); I++){ + N_temp(i,s) = (I->second)/(norm(I->second)); + i++; + } + N = N_temp; + } + + bdry_values normals::operator()(size_t k) const + { + slice s(1,N.rows()); + return bdry_values(N(s,k),rbfs_hash,n,m); + } + + bool normals::operator==(const normals& M) const + { + return rbfs_hash == M.rbfs_hash; + } + + bool normals::operator!=(const normals& M) const + { + return rbfs_hash != M.rbfs_hash; + } + + //*************** 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,v.m); + } + + 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,v.m); + } +}