Mercurial > hg > kwantix
changeset 40:eaa99e09607d
Remove indentation due to namespace scope
| author | Jordi Gutiérrez Hermoso <jordigh@gmail.com> |
|---|---|
| date | Fri, 12 Mar 2010 09:21:07 -0600 |
| parents | 664d2f784d14 |
| children | 834f8e778a65 |
| files | src/bvp.cpp src/ddm.cpp src/diff_op.cpp src/error.cpp src/func.cpp src/include/bvp.hpp src/include/ddm.hpp src/include/diff_op.hpp src/include/error.hpp src/include/func.hpp src/include/interp_values.hpp src/include/interpolator.hpp src/include/linalg.hpp src/include/rbf.hpp src/include/utils.hpp src/include/vtkplot.hpp src/interp_values.cpp src/interpolator.cpp src/linalg.cpp src/main-Laplace.cpp src/main-profiling.cpp src/main-sw-euler.cpp src/main-sw-rk4.cpp src/rbf.cpp src/utils.cpp src/vtkplot.cpp |
| diffstat | 26 files changed, 4742 insertions(+), 4744 deletions(-) [+] |
line wrap: on
line diff
--- a/src/bvp.cpp +++ b/src/bvp.cpp @@ -11,136 +11,121 @@ namespace kwantix { - using std::pair; - using std::make_pair; - - //******************** domain stuff ************************************ - - domain::domain(size_t dimension) - { - dim = dimension; - } - domain::domain(size_t dimension, set<point> intr, - set<point> bdry, map<point, vector> ns) - { - dim = dimension; - add_to_interior(intr); - add_to_boundary(bdry); - add_to_normals(ns); - } - - domain::domain(string intr, string bdry, string ns) - { - - bool intr_empty, bdry_empty; - - matrix intr_m ; - try - { - intr_m = read_matrix(intr); - intr_empty = false; - } - catch(endOfFile) - { - intr_empty = true; - } - - matrix bdry_m; - try - { - bdry_m = read_matrix(bdry); - bdry_empty = false; - } - catch(endOfFile) - { - bdry_empty = true; - } +using std::pair; +using std::make_pair; - matrix ns_m; - try - { - ns_m = read_matrix(ns); - } - catch(endOfFile& exc) - { - if(!bdry_empty) - { - exc.reason = "Boundary is not empty, so normals cannot be either."; - throw exc; - } - } +//******************** domain stuff ************************************ + +domain::domain(size_t dimension) +{ + dim = dimension; +} +domain::domain(size_t dimension, set<point> intr, + set<point> bdry, map<point, vector> ns) +{ + dim = dimension; + add_to_interior(intr); + add_to_boundary(bdry); + add_to_normals(ns); +} + +domain::domain(string intr, string bdry, string ns) +{ + + bool intr_empty, bdry_empty; - dim = intr_m.cols(); - if( bdry_m.cols() != dim and !bdry_empty ) + matrix intr_m ; + try + { + intr_m = read_matrix(intr); + intr_empty = false; + } + catch(endOfFile) + { + intr_empty = true; + } + + matrix bdry_m; + try + { + bdry_m = read_matrix(bdry); + bdry_empty = false; + } + catch(endOfFile) + { + bdry_empty = true; + } + + matrix ns_m; + try + { + ns_m = read_matrix(ns); + } + catch(endOfFile& exc) + { + if(!bdry_empty) { - badArgument exc; - exc.reason = - "Wrong parameters for domain from filename. \n" - "Dimension of boundary (columns) must equal dimension of interior."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - if( ns_m.cols() != 2*dim and !bdry_empty) - { - badArgument exc; - exc.reason = - "Wrong parameters for domain from filename. \n" - "Dimension of normals (columns) must equal twice the \n" - "dimension of the interior and the boundary."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - slice s(1,dim), s1(1,dim), s2(dim+1,2*dim); - - if(!intr_empty) - for(size_t i = 1; i <= intr_m.rows(); i++) - add_to_interior( intr_m(i,s) ); - - if(!bdry_empty){ - for(size_t i = 1; i <= bdry_m.rows(); i++) - add_to_boundary( bdry_m(i,s) ); - for(size_t i = 1; i <= ns_m.rows(); i++) - add_to_normals(ns_m(i,s1), ns_m(i, s2)); + exc.reason = "Boundary is not empty, so normals cannot be either."; + throw exc; } } - domain::~domain() - { - //Nothing! - } - - //This clears any data already in the domain and sets the dimension. - void domain::set_dimension(size_t dimension) + dim = intr_m.cols(); + if( bdry_m.cols() != dim and !bdry_empty ) { - dim = dimension; - interior.clear(); - boundary.clear(); - normals.clear(); + badArgument exc; + exc.reason = + "Wrong parameters for domain from filename. \n" + "Dimension of boundary (columns) must equal dimension of interior."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } - - //Add information to the domain - void domain::add_to_interior(const set<point> &intr) + if( ns_m.cols() != 2*dim and !bdry_empty) { - for(auto I = intr.begin(); I != intr.end(); I++) - { - if(I -> size() != dim) - { - badArgument exc; - exc.reason = - "Cannot assign to domain's interior: inconformant dimensions."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - interior.insert(*I); - } + badArgument exc; + exc.reason = + "Wrong parameters for domain from filename. \n" + "Dimension of normals (columns) must equal twice the \n" + "dimension of the interior and the boundary."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } + slice s(1,dim), s1(1,dim), s2(dim+1,2*dim); + + if(!intr_empty) + for(size_t i = 1; i <= intr_m.rows(); i++) + add_to_interior( intr_m(i,s) ); + + if(!bdry_empty){ + for(size_t i = 1; i <= bdry_m.rows(); i++) + add_to_boundary( bdry_m(i,s) ); + for(size_t i = 1; i <= ns_m.rows(); i++) + add_to_normals(ns_m(i,s1), ns_m(i, s2)); } +} - void domain::add_to_interior(const point &intr) +domain::~domain() +{ + //Nothing! +} + +//This clears any data already in the domain and sets the dimension. +void domain::set_dimension(size_t dimension) +{ + dim = dimension; + interior.clear(); + boundary.clear(); + normals.clear(); +} + +//Add information to the domain +void domain::add_to_interior(const set<point> &intr) +{ + for(auto I = intr.begin(); I != intr.end(); I++) { - if(intr.size() != dim) + if(I -> size() != dim) { badArgument exc; exc.reason = @@ -149,29 +134,29 @@ exc.file = __FILE__; throw exc; } - interior.insert(intr); + interior.insert(*I); } - - void domain::add_to_boundary(const set<point> &bdry) +} + +void domain::add_to_interior(const point &intr) +{ + if(intr.size() != dim) { - for(auto I = bdry.begin(); I != bdry.end(); I++) - { - if(I -> size() != dim) - { - badArgument exc; - exc.reason = - "Cannot assign to domain's boundary: inconformant dimensions."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - boundary.insert(*I); - } + badArgument exc; + exc.reason = + "Cannot assign to domain's interior: inconformant dimensions."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } - - void domain::add_to_boundary(const point &bdry) + interior.insert(intr); +} + +void domain::add_to_boundary(const set<point> &bdry) +{ + for(auto I = bdry.begin(); I != bdry.end(); I++) { - if(bdry.size() != dim) + if(I -> size() != dim) { badArgument exc; exc.reason = @@ -180,36 +165,29 @@ exc.file = __FILE__; throw exc; } - boundary.insert(bdry); + boundary.insert(*I); } +} - void domain::add_to_normals(const map<point, vector> &ns) +void domain::add_to_boundary(const point &bdry) +{ + if(bdry.size() != dim) { - for(auto I = ns.begin(); I != ns.end(); I++) - { - if (!kwantix::contains(boundary, I->first)) - { - badArgument exc; - exc.reason = "Bad normal given: must match a point on the boundary."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - if(I->first.size() != dim or I->second.size() != dim) - { - badArgument exc; - exc.reason = "Bad normal given: inconformant dimensions."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - normals.insert(*I); - } + badArgument exc; + exc.reason = + "Cannot assign to domain's boundary: inconformant dimensions."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } + boundary.insert(bdry); +} - void domain::add_to_normals(const point &bdry, const vector &n) +void domain::add_to_normals(const map<point, vector> &ns) +{ + for(auto I = ns.begin(); I != ns.end(); I++) { - if (!kwantix::contains(boundary, bdry)) + if (!kwantix::contains(boundary, I->first)) { badArgument exc; exc.reason = "Bad normal given: must match a point on the boundary."; @@ -217,7 +195,7 @@ exc.file = __FILE__; throw exc; } - if(bdry.size() != dim or n.size() != dim) + if(I->first.size() != dim or I->second.size() != dim) { badArgument exc; exc.reason = "Bad normal given: inconformant dimensions."; @@ -225,118 +203,140 @@ exc.file = __FILE__; throw exc; } - pair<point, vector> pn = make_pair(bdry,n); - normals.insert(pn); - } - - //Read that information - size_t domain::get_dimension() const - { - return dim; - } - const set<point>& domain::get_interior() const - { - return interior; - } - const set<point>& domain::get_boundary() const - { - return boundary; + normals.insert(*I); } - const map<point, vector>& domain::get_normals() const - { - return normals; - } - - //Is point in this domain, whether interior or boundary? - bool domain::contains(const point& p) const +} + +void domain::add_to_normals(const point &bdry, const vector &n) +{ + if (!kwantix::contains(boundary, bdry)) { - if(kwantix::contains(interior, p)) - return true; - if(kwantix::contains(boundary, p)) - return true; - return false; + badArgument exc; + exc.reason = "Bad normal given: must match a point on the boundary."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } - - //**************** BVP stuff ********************************* - - shared_ptr<const domain> BVP::get_domain() const + if(bdry.size() != dim or n.size() != dim) { - return Omega; - } - shared_ptr<const diff_op> BVP::get_diff_op() const - { - return L; - } - shared_ptr<const bdry_diff_op> BVP::get_bdry_diff_op() const - { - return B; - } - const map<point, double>& BVP::get_f() const - { - return f; - } - const map<point, double>& BVP::get_g() const - { - return g; + badArgument exc; + exc.reason = "Bad normal given: inconformant dimensions."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } + pair<point, vector> pn = make_pair(bdry,n); + normals.insert(pn); +} - void BVP::set_f(const map<point, double>& f_in) +//Read that information +size_t domain::get_dimension() const +{ + return dim; +} +const set<point>& domain::get_interior() const +{ + return interior; +} +const set<point>& domain::get_boundary() const +{ + return boundary; +} +const map<point, vector>& domain::get_normals() const +{ + return normals; +} + +//Is point in this domain, whether interior or boundary? +bool domain::contains(const point& p) const +{ + if(kwantix::contains(interior, p)) + return true; + if(kwantix::contains(boundary, p)) + return true; + return false; +} + +//**************** BVP stuff ********************************* + +shared_ptr<const domain> BVP::get_domain() const +{ + return Omega; +} +shared_ptr<const diff_op> BVP::get_diff_op() const +{ + return L; +} +shared_ptr<const bdry_diff_op> BVP::get_bdry_diff_op() const +{ + return B; +} +const map<point, double>& BVP::get_f() const +{ + return f; +} +const map<point, double>& BVP::get_g() const +{ + return g; +} + +void BVP::set_f(const map<point, double>& f_in) +{ + for( auto I = f_in.begin(); I != f_in.end(); I++) { - for( auto I = f_in.begin(); I != f_in.end(); I++) + if( !contains(Omega->get_interior(), I->first ) ) { - if( !contains(Omega->get_interior(), I->first ) ) - { - badArgument exc; - exc.reason = - "Bad specification of f in BVP: " - "gave a point not in the interior."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - f[I->first] = I->second; + badArgument exc; + exc.reason = + "Bad specification of f in BVP: " + "gave a point not in the interior."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } + f[I->first] = I->second; } +} - void BVP::set_g(const map<point, double>& g_in) - { - for( auto I = g_in.begin(); I != g_in.end(); I++) - { - if( !contains(Omega->get_boundary(), I->first ) ) - { - badArgument exc; - exc.reason = - "Bad specification of g in BVP: " - "gave a point not on the boundary."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - g[I->first] = I -> second; - } - } - void BVP::set_f(const realfunc &f_in) +void BVP::set_g(const map<point, double>& g_in) +{ + for( auto I = g_in.begin(); I != g_in.end(); I++) { - for(auto I = Omega->get_interior().begin(); I != Omega->get_interior().end(); I++) - { - f[*I] = f_in(*I); - } - } - - void BVP::set_g(const realfunc &g_in) - { - for(auto I = Omega->get_boundary().begin(); I != Omega->get_boundary().end(); I++) + if( !contains(Omega->get_boundary(), I->first ) ) { - g[*I] = g_in(*I); + badArgument exc; + exc.reason = + "Bad specification of g in BVP: " + "gave a point not on the boundary."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } + g[I->first] = I -> second; } +} +void BVP::set_f(const realfunc &f_in) +{ + for(auto I = Omega->get_interior().begin(); I != Omega->get_interior().end(); I++) + { + f[*I] = f_in(*I); + } +} - shared_ptr<const linear_diff_op2> linear_BVP2::get_linear_diff_op2() const +void BVP::set_g(const realfunc &g_in) +{ + for(auto I = Omega->get_boundary().begin(); I != Omega->get_boundary().end(); I++) { - return boost:: - dynamic_pointer_cast<const linear_diff_op2>(BVP::get_diff_op()); + g[*I] = g_in(*I); } +} + +shared_ptr<const linear_diff_op2> linear_BVP2::get_linear_diff_op2() const +{ + return boost:: + dynamic_pointer_cast<const linear_diff_op2>(BVP::get_diff_op()); +} }
--- a/src/ddm.cpp +++ b/src/ddm.cpp @@ -14,430 +14,430 @@ namespace kwantix { - using namespace std; - using kwantix::vector; - using boost::shared_ptr; - using boost::dynamic_pointer_cast; +using namespace std; +using kwantix::vector; +using boost::shared_ptr; +using boost::dynamic_pointer_cast; - //************* ddm_bdry_diff_op stuff ***************** +//************* ddm_bdry_diff_op stuff ***************** - ddm_bdry_diff_op::ddm_bdry_diff_op(shared_ptr<const bdry_diff_op> B_in, - shared_ptr<const bdry_diff_op> Bprime_in, - const set<point>& ibps) - { - B = B_in ; - intr_bdry_pts = ibps; - Bprime = Bprime_in; - } +ddm_bdry_diff_op::ddm_bdry_diff_op(shared_ptr<const bdry_diff_op> B_in, + shared_ptr<const bdry_diff_op> Bprime_in, + const set<point>& ibps) +{ + B = B_in ; + intr_bdry_pts = ibps; + Bprime = Bprime_in; +} - double ddm_bdry_diff_op::at(const realfunc &f, const point &p) const - { - if( kwantix::contains(intr_bdry_pts, p)) - return Bprime -> at(f,p); +double ddm_bdry_diff_op::at(const realfunc &f, const point &p) const +{ + if( kwantix::contains(intr_bdry_pts, p)) + return Bprime -> at(f,p); - return B -> at(f,p); - } + return B -> at(f,p); +} + +double ddm_bdry_diff_op::at(const realfunc &f, const point &p, + const vector &n) const +{ + if( kwantix::contains(intr_bdry_pts, p)) + return Bprime -> at(f,p,n); + + return B -> at(f,p,n); +} - double ddm_bdry_diff_op::at(const realfunc &f, const point &p, - const vector &n) const +// ************ ddm stuff ****************************** + +ddm::ddm(const set<shared_ptr<const domain> >& ds, + shared_ptr<const BVP> thebvp) +{ + //Gotta check this is actually a domain decomposition... + set<point> union_interior; + set<point> union_boundary; + + for(auto i = ds.begin(); i != ds.end(); i++) { - if( kwantix::contains(intr_bdry_pts, p)) - return Bprime -> at(f,p,n); - - return B -> at(f,p,n); + set<point> intr = (*i) -> get_interior(); + set<point> bdry = (*i) -> get_boundary(); + union_interior.insert(intr.begin(), intr.end()); + union_boundary.insert(bdry.begin(), bdry.end()); + } + + bvp = thebvp; + shared_ptr<const domain> Omega = bvp -> get_domain(); + set<point> interior = Omega -> get_interior(); + set<point> boundary = Omega -> get_boundary(); + + if( interior != union_interior) + { + badArgument exc; + exc.reason = + "Bad argument in domain decomposition method constructor: \n" + "The union of the interior of the proposed domains does not \n" + "equal the interior of the domain."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } + + if(!kwantix::includes(boundary,union_boundary) ) + { + badArgument exc; + exc.reason = + "Bad argument in domain decomposition method constructor: \n" + "The union of the boundary of the proposed domains does not \n" + "contain the boundary of the domain."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } - - // ************ ddm stuff ****************************** + + domains = ds; + tolerance = 1e-6; + +} + +ddm::~ddm() +{ + //Nothing! +} + +void ddm::set_tolerance(double tol) +{ + tolerance = tol; +} + +double ddm::operator()(const point& p) const +{ + return at(p); +} - ddm::ddm(const set<shared_ptr<const domain> >& ds, - shared_ptr<const BVP> thebvp) +// ************** additive_schwarz_ddm stuff ************************* + +template <typename RBF> additive_schwarz_ddm<RBF>:: +additive_schwarz_ddm(const set<shared_ptr<const domain> >& ds, + const shared_ptr<const linear_BVP2> thebvp) : + ddm(ds,thebvp) +{ + + shared_ptr<const linear_diff_op2> + L = dynamic_pointer_cast< + const linear_diff_op2>(this -> bvp -> get_diff_op()); + + shared_ptr<const bdry_diff_op> + B = this -> bvp -> get_bdry_diff_op(); + + map<point, double> global_f = this -> bvp -> get_f(); + map<point, double> global_g = this -> bvp -> get_g(); + + //Define a bvp for each domain and assign it an interpolator. + for(auto i = this -> domains.begin(); i != this -> domains.end(); i++) { - //Gotta check this is actually a domain decomposition... - set<point> union_interior; - set<point> union_boundary; + shared_ptr<const overlapping_domain> this_domain = + dynamic_pointer_cast<const overlapping_domain>(*i); + + //Define the interiors and f's + set<point> this_intr = this_domain -> get_interior(); + map<point, double> this_f; + for(auto spi = this_intr.begin(); spi != this_intr.end(); spi++) + this_f[*spi] = global_f[*spi]; - for(auto i = ds.begin(); i != ds.end(); i++) + //Define the boundaries and g's + set<point> this_bdry = this_domain -> get_boundary(); + map<point, double> this_g; + set<point> interior_boundary_pts; + for(auto spi = this_bdry.begin(); spi != this_bdry.end(); spi++) { - set<point> intr = (*i) -> get_interior(); - set<point> bdry = (*i) -> get_boundary(); - union_interior.insert(intr.begin(), intr.end()); - union_boundary.insert(bdry.begin(), bdry.end()); + if(this_domain -> which_domain(*spi).get() == 0) + { + this_g[*spi] = global_g[*spi]; + } + else + { + interior_boundary_pts.insert(*spi); + this_g[*spi] = 0; //Init to zero artificial boundary conditions. + } } - bvp = thebvp; - shared_ptr<const domain> Omega = bvp -> get_domain(); - set<point> interior = Omega -> get_interior(); - set<point> boundary = Omega -> get_boundary(); + //Define the boundary operator, using Dirichlet for artificial + //boundaries + shared_ptr<dirichlet_op> D(new dirichlet_op); + shared_ptr<ddm_bdry_diff_op> + this_B( new ddm_bdry_diff_op(B, D, interior_boundary_pts) ); + + //Put all this info into an interpolator. + shared_ptr<linear_BVP2> + this_bvp(new linear_BVP2(this_domain, L, this_B, this_f, this_g) ); + shared_ptr<interpolator<RBF> > + rbf_ptr( new interpolator<RBF>(this_bvp) ); + phis[this_domain] = rbf_ptr; + } + + //Apply actual additive Schwarz DDM algorithm here + solve(); +} + +template <typename RBF> +double additive_schwarz_ddm<RBF>::at(const point& p) const +{ + //If p is one of the domain points, then it's used. If more than + //one domain contains p, then the average is taken. If no domain + //contains p, we exhaustively search for the closest point in the + //domains and use the interpolator from that domain to evaluate. + + set<shared_ptr<const interpolator<RBF> > > + relevant_interpolators = which_interps(p); + + if(relevant_interpolators.size() != 0)//p is in one of the + //domains. + return avg_interp(relevant_interpolators, p); + + + //Else, p is not one of the grid points. Find closest grid point + //and evaluate on the domain(s) that the point belongs to. + + + //Uh, just begin with closest being some point in the whatever + //domain. + point c = *(((*( this->domains.begin() )) + -> get_interior()).begin() ); + + //Search each domain's interior and boundary. Can't do better than + //exhaustive search. + for(auto i = this -> domains.begin(); i != this -> domains.end(); i++) + { + for(set<point>::iterator j = (*i) -> get_interior().begin(); + j != (*i) -> get_interior().end(); j++) + if( norm(*j - p) < norm(c - p)) + c = *j; + + for(set<point>::iterator j = (*i) -> get_boundary().begin(); + j != (*i) -> get_boundary().end(); j++) + if( norm(*j - p) < norm(c - p)) + c = *j; + } + + return avg_interp( which_interps(c), p); + +} + +template <typename RBF> +set<shared_ptr<const interpolator<RBF> > > additive_schwarz_ddm<RBF>:: +which_interps(const point& p) const +{ + set<shared_ptr<const interpolator<RBF> > > relevant_interpolators; + for(set<shared_ptr<const domain> >::iterator i = this -> domains.begin(); + i != this -> domains.end(); i++){ + if( (*i) -> contains(p)){ + shared_ptr<const overlapping_domain> j = + dynamic_pointer_cast<const overlapping_domain>(*i); + relevant_interpolators.insert(phis.at(j)); + //at is not in current STL standard; but it is necessary here + //because operator[] can't be used here since this is a const + //function. at is currently a GNU extension, and we're using + //it. + } + } + return relevant_interpolators; +} + +template <typename RBF> +double additive_schwarz_ddm<RBF>:: +avg_interp(set<shared_ptr<const interpolator<RBF> > > relevant_interpolators, + const point& p) const +{ + double result = 0; + int n = 0; + for(typename set<shared_ptr<const interpolator<RBF> > >:: + iterator i = relevant_interpolators.begin() ; + i != relevant_interpolators.end(); i++){ + result += (*i) -> at(p); + n++; + } + return result/n; +} - if( interior != union_interior) +template <typename RBF> +void additive_schwarz_ddm<RBF>::solve() +{ + //Recap: each domain already has an interpolator associated to it + //and all overlapping domain information has been defined. Just + //need to iterate on the boundary conditions. Method converges + //when 2-norm (Frobenius, if it were a matrix) of the artificial + //boundary does not change by more than tolerance. + using std::make_pair; + + vector newv = at_all_points(); + vector oldv(newv.size()); + double change; + do{ + oldv = newv; + + //Each domain will have new g's. + map<shared_ptr<const overlapping_domain>, map<point, double> > + new_bdry_assts; + for(set<shared_ptr<const domain> >::iterator + i = this -> domains.begin(); i != this -> domains.end(); i++){ + shared_ptr<const overlapping_domain> d = + dynamic_pointer_cast<const overlapping_domain>(*i); + for(set<point>::iterator j = (*i) -> get_boundary().begin(); + j != (*i) -> get_boundary().end(); j++) + if( d -> which_domain(*j).get() != 0){ + + + new_bdry_assts[d]. + insert(make_pair(*j, + phis[d -> which_domain(*j)] -> at(*j)) + ); + } + } + + //Now assign to each interpolator the modified boundary. + for(typename map<shared_ptr<const overlapping_domain>, + shared_ptr<interpolator<RBF> > >:: + iterator i = phis.begin(); i != phis.end(); i++) + i -> second -> set_g(new_bdry_assts[i -> first]); + + newv = at_all_points(); + change = norm(oldv-newv); + + //debug + cout << "Change: " << change << endl; + }while( change > this -> tolerance); +} + +//Evaluate problem at all artificial boundary points. +template <typename RBF> +vector additive_schwarz_ddm<RBF>::at_all_points() const +{ + set<point> art_bdry; + for(set<shared_ptr<const domain> >::const_iterator i = + this -> domains.begin(); i != this -> domains.end(); i++) + { + shared_ptr<const overlapping_domain> j = + dynamic_pointer_cast<const overlapping_domain>(*i), + zero; //The zero pointer. + + for(set<point>::iterator p = j -> get_boundary().begin(); + p != j -> get_boundary().end(); p++) { + if(j -> which_domain(*p) != zero) + art_bdry.insert(*p); + } + } + + set<point>::iterator I = art_bdry.begin(); + vector result(art_bdry.size() ); + for(size_t i = 1; i <= result.size(); i++){ + result(i) = at(*I); + I++; + } + return result; +} + +//************** overlapping_domain stuff ************* + + +overlapping_domain:: +overlapping_domain(string intr, string bdry, string ns) + :domain(intr,bdry,ns) +{ +} + +overlapping_domain:: +overlapping_domain(string intr, string bdry, string ns, + const set<shared_ptr<const overlapping_domain> >& ols, + const map<point, + shared_ptr<const overlapping_domain> >& bdry_asst) + :domain(intr,bdry,ns) +{ + set<point> bdry_copy = get_boundary(); + for(map<point, shared_ptr<const overlapping_domain> >::const_iterator + i = bdry_asst.begin(); i != bdry_asst.end(); i++) + { + if(!kwantix::contains(ols, i->second) or + !kwantix::contains(bdry_copy, i->first)){ badArgument exc; exc.reason = - "Bad argument in domain decomposition method constructor: \n" - "The union of the interior of the proposed domains does not \n" - "equal the interior of the domain."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - - if(!kwantix::includes(boundary,union_boundary) ) - { - badArgument exc; - exc.reason = - "Bad argument in domain decomposition method constructor: \n" - "The union of the boundary of the proposed domains does not \n" - "contain the boundary of the domain."; + "Bad argument in overlapping_domain constructor: \n" + "Every boundary assignment must be a boundary point \n" + "assigned to some overlapping domain."; exc.line = __LINE__; exc.file = __FILE__; throw exc; } - - domains = ds; - tolerance = 1e-6; - } - - ddm::~ddm() - { - //Nothing! - } - - void ddm::set_tolerance(double tol) - { - tolerance = tol; - } - - double ddm::operator()(const point& p) const - { - return at(p); - } - - // ************** additive_schwarz_ddm stuff ************************* - - template <typename RBF> additive_schwarz_ddm<RBF>:: - additive_schwarz_ddm(const set<shared_ptr<const domain> >& ds, - const shared_ptr<const linear_BVP2> thebvp) : - ddm(ds,thebvp) - { - - shared_ptr<const linear_diff_op2> - L = dynamic_pointer_cast< - const linear_diff_op2>(this -> bvp -> get_diff_op()); + overlappers = ols; + boundary_assignments = bdry_asst; +} - shared_ptr<const bdry_diff_op> - B = this -> bvp -> get_bdry_diff_op(); - - map<point, double> global_f = this -> bvp -> get_f(); - map<point, double> global_g = this -> bvp -> get_g(); - - //Define a bvp for each domain and assign it an interpolator. - for(auto i = this -> domains.begin(); i != this -> domains.end(); i++) - { - shared_ptr<const overlapping_domain> this_domain = - dynamic_pointer_cast<const overlapping_domain>(*i); - - //Define the interiors and f's - set<point> this_intr = this_domain -> get_interior(); - map<point, double> this_f; - for(auto spi = this_intr.begin(); spi != this_intr.end(); spi++) - this_f[*spi] = global_f[*spi]; - - //Define the boundaries and g's - set<point> this_bdry = this_domain -> get_boundary(); - map<point, double> this_g; - set<point> interior_boundary_pts; - for(auto spi = this_bdry.begin(); spi != this_bdry.end(); spi++) - { - if(this_domain -> which_domain(*spi).get() == 0) - { - this_g[*spi] = global_g[*spi]; - } - else - { - interior_boundary_pts.insert(*spi); - this_g[*spi] = 0; //Init to zero artificial boundary conditions. - } - } - - //Define the boundary operator, using Dirichlet for artificial - //boundaries - shared_ptr<dirichlet_op> D(new dirichlet_op); - shared_ptr<ddm_bdry_diff_op> - this_B( new ddm_bdry_diff_op(B, D, interior_boundary_pts) ); - - //Put all this info into an interpolator. - shared_ptr<linear_BVP2> - this_bvp(new linear_BVP2(this_domain, L, this_B, this_f, this_g) ); - shared_ptr<interpolator<RBF> > - rbf_ptr( new interpolator<RBF>(this_bvp) ); - phis[this_domain] = rbf_ptr; - } - - //Apply actual additive Schwarz DDM algorithm here - solve(); - } - - template <typename RBF> - double additive_schwarz_ddm<RBF>::at(const point& p) const - { - //If p is one of the domain points, then it's used. If more than - //one domain contains p, then the average is taken. If no domain - //contains p, we exhaustively search for the closest point in the - //domains and use the interpolator from that domain to evaluate. - - set<shared_ptr<const interpolator<RBF> > > - relevant_interpolators = which_interps(p); - - if(relevant_interpolators.size() != 0)//p is in one of the - //domains. - return avg_interp(relevant_interpolators, p); +overlapping_domain::overlapping_domain(size_t dimension) + :domain(dimension){} +overlapping_domain:: +overlapping_domain(size_t dimension, set<point> intr, + set<point> bdry, map<point, vector> ns) + :domain(dimension, intr, bdry, ns){} - //Else, p is not one of the grid points. Find closest grid point - //and evaluate on the domain(s) that the point belongs to. - - - //Uh, just begin with closest being some point in the whatever - //domain. - point c = *(((*( this->domains.begin() )) - -> get_interior()).begin() ); - - //Search each domain's interior and boundary. Can't do better than - //exhaustive search. - for(auto i = this -> domains.begin(); i != this -> domains.end(); i++) - { - for(set<point>::iterator j = (*i) -> get_interior().begin(); - j != (*i) -> get_interior().end(); j++) - if( norm(*j - p) < norm(c - p)) - c = *j; - - for(set<point>::iterator j = (*i) -> get_boundary().begin(); - j != (*i) -> get_boundary().end(); j++) - if( norm(*j - p) < norm(c - p)) - c = *j; - } - - return avg_interp( which_interps(c), p); - - } - - template <typename RBF> - set<shared_ptr<const interpolator<RBF> > > additive_schwarz_ddm<RBF>:: - which_interps(const point& p) const - { - set<shared_ptr<const interpolator<RBF> > > relevant_interpolators; - for(set<shared_ptr<const domain> >::iterator i = this -> domains.begin(); - i != this -> domains.end(); i++){ - if( (*i) -> contains(p)){ - shared_ptr<const overlapping_domain> j = - dynamic_pointer_cast<const overlapping_domain>(*i); - relevant_interpolators.insert(phis.at(j)); - //at is not in current STL standard; but it is necessary here - //because operator[] can't be used here since this is a const - //function. at is currently a GNU extension, and we're using - //it. - } - } - return relevant_interpolators; - } - - template <typename RBF> - double additive_schwarz_ddm<RBF>:: - avg_interp(set<shared_ptr<const interpolator<RBF> > > relevant_interpolators, - const point& p) const - { - double result = 0; - int n = 0; - for(typename set<shared_ptr<const interpolator<RBF> > >:: - iterator i = relevant_interpolators.begin() ; - i != relevant_interpolators.end(); i++){ - result += (*i) -> at(p); - n++; - } - return result/n; - } - - template <typename RBF> - void additive_schwarz_ddm<RBF>::solve() - { - //Recap: each domain already has an interpolator associated to it - //and all overlapping domain information has been defined. Just - //need to iterate on the boundary conditions. Method converges - //when 2-norm (Frobenius, if it were a matrix) of the artificial - //boundary does not change by more than tolerance. - using std::make_pair; - - vector newv = at_all_points(); - vector oldv(newv.size()); - double change; - do{ - oldv = newv; - - //Each domain will have new g's. - map<shared_ptr<const overlapping_domain>, map<point, double> > - new_bdry_assts; - for(set<shared_ptr<const domain> >::iterator - i = this -> domains.begin(); i != this -> domains.end(); i++){ - shared_ptr<const overlapping_domain> d = - dynamic_pointer_cast<const overlapping_domain>(*i); - for(set<point>::iterator j = (*i) -> get_boundary().begin(); - j != (*i) -> get_boundary().end(); j++) - if( d -> which_domain(*j).get() != 0){ - - - new_bdry_assts[d]. - insert(make_pair(*j, - phis[d -> which_domain(*j)] -> at(*j)) - ); - } - } - - //Now assign to each interpolator the modified boundary. - for(typename map<shared_ptr<const overlapping_domain>, - shared_ptr<interpolator<RBF> > >:: - iterator i = phis.begin(); i != phis.end(); i++) - i -> second -> set_g(new_bdry_assts[i -> first]); - - newv = at_all_points(); - change = norm(oldv-newv); - - //debug - cout << "Change: " << change << endl; - }while( change > this -> tolerance); - } +set<shared_ptr<const overlapping_domain> > +overlapping_domain::get_domains() const +{ + return overlappers; +} - //Evaluate problem at all artificial boundary points. - template <typename RBF> - vector additive_schwarz_ddm<RBF>::at_all_points() const +shared_ptr<const overlapping_domain> +overlapping_domain::which_domain(const point& p) const +{ + if( !kwantix::contains(boundary_assignments, p) ) { - set<point> art_bdry; - for(set<shared_ptr<const domain> >::const_iterator i = - this -> domains.begin(); i != this -> domains.end(); i++) - { - shared_ptr<const overlapping_domain> j = - dynamic_pointer_cast<const overlapping_domain>(*i), - zero; //The zero pointer. - - for(set<point>::iterator p = j -> get_boundary().begin(); - p != j -> get_boundary().end(); p++) - { - if(j -> which_domain(*p) != zero) - art_bdry.insert(*p); - } - } + shared_ptr<const overlapping_domain> zero; + return zero; + } + return boundary_assignments.at(p); +} + +void +overlapping_domain:: +set_overlapper_info(const point& p, const shared_ptr<overlapping_domain> o) +{ + if(kwantix::contains(this -> get_boundary(), p)) + boundary_assignments[p] = o; +} - set<point>::iterator I = art_bdry.begin(); - vector result(art_bdry.size() ); - for(size_t i = 1; i <= result.size(); i++){ - result(i) = at(*I); - I++; - } - return result; - } - - //************** overlapping_domain stuff ************* - +//************** friends of overlapping_domain ******************* - overlapping_domain:: - overlapping_domain(string intr, string bdry, string ns) - :domain(intr,bdry,ns) +void +set_overlapper_info(set<shared_ptr<overlapping_domain> > domains) +{ + for(auto d = domains.begin(); d != domains.end(); d++) { - } - - overlapping_domain:: - overlapping_domain(string intr, string bdry, string ns, - const set<shared_ptr<const overlapping_domain> >& ols, - const map<point, - shared_ptr<const overlapping_domain> >& bdry_asst) - :domain(intr,bdry,ns) - { - set<point> bdry_copy = get_boundary(); - for(map<point, shared_ptr<const overlapping_domain> >::const_iterator - i = bdry_asst.begin(); i != bdry_asst.end(); i++) + for(auto p = (*d) -> get_boundary().begin(); + p != (*d) -> get_boundary().end(); p++) { - if(!kwantix::contains(ols, i->second) or - !kwantix::contains(bdry_copy, i->first)){ - badArgument exc; - exc.reason = - "Bad argument in overlapping_domain constructor: \n" - "Every boundary assignment must be a boundary point \n" - "assigned to some overlapping domain."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - } - overlappers = ols; - boundary_assignments = bdry_asst; - } - - overlapping_domain::overlapping_domain(size_t dimension) - :domain(dimension){} - overlapping_domain:: - overlapping_domain(size_t dimension, set<point> intr, - set<point> bdry, map<point, vector> ns) - :domain(dimension, intr, bdry, ns){} - - - set<shared_ptr<const overlapping_domain> > - overlapping_domain::get_domains() const - { - return overlappers; - } - - shared_ptr<const overlapping_domain> - overlapping_domain::which_domain(const point& p) const - { - if( !kwantix::contains(boundary_assignments, p) ) - { - shared_ptr<const overlapping_domain> zero; - return zero; - } - return boundary_assignments.at(p); - } - - void - overlapping_domain:: - set_overlapper_info(const point& p, const shared_ptr<overlapping_domain> o) - { - if(kwantix::contains(this -> get_boundary(), p)) - boundary_assignments[p] = o; - } - - //************** friends of overlapping_domain ******************* - - void - set_overlapper_info(set<shared_ptr<overlapping_domain> > domains) - { - for(auto d = domains.begin(); d != domains.end(); d++) - { - for(auto p = (*d) -> get_boundary().begin(); - p != (*d) -> get_boundary().end(); p++) + for(auto d_other = domains.begin(); d_other != domains.end(); d_other++) { - for(auto d_other = domains.begin(); d_other != domains.end(); d_other++) - { - if( - kwantix::contains((*d_other ) -> get_interior(), *p) - ){ - (*d) -> boundary_assignments[*p] = *d_other; - (*d) -> overlappers.insert(*d_other); - break; //FIXME: We're assuming no three domains overlap at - //one point. - } + if( + kwantix::contains((*d_other ) -> get_interior(), *p) + ){ + (*d) -> boundary_assignments[*p] = *d_other; + (*d) -> overlappers.insert(*d_other); + break; //FIXME: We're assuming no three domains overlap at + //one point. } } } } +} - //Instantiations - template class additive_schwarz_ddm<piecewise_polynomial>; - template class additive_schwarz_ddm<thin_plate_spline>; - template class additive_schwarz_ddm<multiquadric>; - template class additive_schwarz_ddm<inverse_multiquadric>; - template class additive_schwarz_ddm<inverse_quadratic>; - template class additive_schwarz_ddm<gaussian>; +//Instantiations +template class additive_schwarz_ddm<piecewise_polynomial>; +template class additive_schwarz_ddm<thin_plate_spline>; +template class additive_schwarz_ddm<multiquadric>; +template class additive_schwarz_ddm<inverse_multiquadric>; +template class additive_schwarz_ddm<inverse_quadratic>; +template class additive_schwarz_ddm<gaussian>; }
--- a/src/diff_op.cpp +++ b/src/diff_op.cpp @@ -4,89 +4,89 @@ namespace kwantix{ - using namespace std; +using namespace std; - //************** Differential operator functions ********************* - double diff_op::operator()(const realfunc &f, const point &p) const - { - return at(f,p); - } +//************** Differential operator functions ********************* +double diff_op::operator()(const realfunc &f, const point &p) const +{ + return at(f,p); +} - double linear_diff_op2::operator()(const realfunc &f, const point - &p) const - { - return diff_op::operator()(f,p); - } +double linear_diff_op2::operator()(const realfunc &f, const point + &p) const +{ + return diff_op::operator()(f,p); +} - double bdry_diff_op::operator()(const realfunc &f, - const point &p) const - { - return at(f,p); - } - double bdry_diff_op::operator()(const realfunc &f, - const point &p, const vector &n) - const - { - return at(f,p,n); - } +double bdry_diff_op::operator()(const realfunc &f, + const point &p) const +{ + return at(f,p); +} +double bdry_diff_op::operator()(const realfunc &f, + const point &p, const vector &n) + const +{ + return at(f,p,n); +} - double dirichlet_op::at(const realfunc &f, - const point &p) const - { - return f(p); - } - double dirichlet_op::at(const realfunc &f, const point &p, - const vector &n) const - { - n.size(); - return f(p); - } +double dirichlet_op::at(const realfunc &f, + const point &p) const +{ + return f(p); +} +double dirichlet_op::at(const realfunc &f, const point &p, + const vector &n) const +{ + n.size(); + return f(p); +} - double neumann_op::at(const realfunc &f, const point &p, - const vector &n) const - { - size_t dim = n.size(); - vector grad(dim); - for(size_t i = 1; i <= dim; i++) - grad(i) =f.d(p,i); - return grad%n/norm(n); - } +double neumann_op::at(const realfunc &f, const point &p, + const vector &n) const +{ + size_t dim = n.size(); + vector grad(dim); + for(size_t i = 1; i <= dim; i++) + grad(i) =f.d(p,i); + return grad%n/norm(n); +} - double Id_op::at(const realfunc &f, const point &p) const - { - return f(p); - } +double Id_op::at(const realfunc &f, const point &p) const +{ + return f(p); +} + +double del1::at(const realfunc &f, const point &p) const +{ + return f.d(p,1); +} + +double del1::at(const realfunc &f, const point &p, size_t i) const +{ + return f.d(p,i); +} - double del1::at(const realfunc &f, const point &p) const - { - return f.d(p,1); - } +double del2::at(const realfunc &f, const point &p) const +{ + return f.d2(p,1,1); +} - double del1::at(const realfunc &f, const point &p, size_t i) const - { - return f.d(p,i); - } +double del2::at(const realfunc &f, const point &p, + size_t i, size_t j) const +{ + return f.d2(p,i,j); +} - double del2::at(const realfunc &f, const point &p) const - { - return f.d2(p,1,1); +double Laplacian::at(const realfunc &f, const point &p) const +{ + size_t n = p.size(); + double result = 0; + del2 D2; + for (size_t i = 1; i <=n; i++){ + result += D2.at(f, p, i, i); } - - double del2::at(const realfunc &f, const point &p, - size_t i, size_t j) const - { - return f.d2(p,i,j); - } - - double Laplacian::at(const realfunc &f, const point &p) const - { - size_t n = p.size(); - double result = 0; - del2 D2; - for (size_t i = 1; i <=n; i++){ - result += D2.at(f, p, i, i); - } - return result; - } + return result; } +}
--- a/src/error.cpp +++ b/src/error.cpp @@ -2,88 +2,88 @@ #include "include/error.hpp" namespace kwantix{ - void errorHandler(const char * reason, const char * file, - int line, int gsl_errno){ - //This exception is so common that we will want more information - //when it happens. - if(reason == string("index out of range")){ //GSL gives this message. - throw indexOutOfRange(reason,file,line); +void errorHandler(const char * reason, const char * file, + int line, int gsl_errno){ + //This exception is so common that we will want more information + //when it happens. + if(reason == string("index out of range")){ //GSL gives this message. + throw indexOutOfRange(reason,file,line); + } + + else{ //Other exceptions are more generic + + switch(gsl_errno){ + case -2: + throw noConvergence(reason,file,line); break; + case 1: + throw badDomain(reason,file,line); break; + case 2: + throw badRange(reason,file,line); break; + case 3: + throw badPointer(reason,file,line); break; + case 4: + throw badArgument(reason,file,line); break; + case 5: + throw failure(reason,file,line); break; + case 6: + throw failedFactorisation(reason,file,line); break; + case 7: + throw failedSanity(reason,file,line); break; + case 8: + throw outOfMemory(reason,file,line); break; + case 9: + throw badFunction(reason,file,line); break; + case 10: + throw runAway(reason,file,line); break; + case 11: + throw maxIterations(reason,file,line); break; + case 12: + throw divideByZero(reason,file,line); break; + case 13: + throw badTolerance(reason,file,line); break; + case 14: + throw aboveTolerance(reason,file,line); break; + case 15: + throw underflow(reason,file,line); break; + case 16: + throw overflow(reason,file,line); break; + case 17: + throw lossOfAccuracy(reason,file,line); break; + case 18: + throw roundOffError(reason,file,line); break; + case 19: + throw inconformantSizes(reason,file,line); break; + case 20: + throw matrixNotSquare(reason,file,line); break; + case 21: + throw singularityFound(reason,file,line); break; + case 22: + throw integralOrSeriesDivergent(reason,file,line); break; + case 23: + throw badHardware(reason,file,line); break; + case 24: + throw notImplemented(reason,file,line); break; + case 25: + throw cacheLimitExceeded(reason,file,line); break; + case 26: + throw tableLimitExceeded(reason,file,line); break; + case 27: + throw iterationNotProgressing(reason,file,line); break; + case 28: + throw jacobiansNotImprovingSolution(reason,file,line); break; + case 29: + cannotReachToleranceInF(reason,file,line); break; + case 30: + throw cannotReachToleranceInX(reason,file,line); break; + case 31: + throw cannotReachToleranceInGradient(reason,file,line); break; + case 32: + throw endOfFile(reason,file,line); break; + default: //Corresponds to GSL_ERRNO=-1. + throw error(reason,file,line); } - else{ //Other exceptions are more generic - - switch(gsl_errno){ - case -2: - throw noConvergence(reason,file,line); break; - case 1: - throw badDomain(reason,file,line); break; - case 2: - throw badRange(reason,file,line); break; - case 3: - throw badPointer(reason,file,line); break; - case 4: - throw badArgument(reason,file,line); break; - case 5: - throw failure(reason,file,line); break; - case 6: - throw failedFactorisation(reason,file,line); break; - case 7: - throw failedSanity(reason,file,line); break; - case 8: - throw outOfMemory(reason,file,line); break; - case 9: - throw badFunction(reason,file,line); break; - case 10: - throw runAway(reason,file,line); break; - case 11: - throw maxIterations(reason,file,line); break; - case 12: - throw divideByZero(reason,file,line); break; - case 13: - throw badTolerance(reason,file,line); break; - case 14: - throw aboveTolerance(reason,file,line); break; - case 15: - throw underflow(reason,file,line); break; - case 16: - throw overflow(reason,file,line); break; - case 17: - throw lossOfAccuracy(reason,file,line); break; - case 18: - throw roundOffError(reason,file,line); break; - case 19: - throw inconformantSizes(reason,file,line); break; - case 20: - throw matrixNotSquare(reason,file,line); break; - case 21: - throw singularityFound(reason,file,line); break; - case 22: - throw integralOrSeriesDivergent(reason,file,line); break; - case 23: - throw badHardware(reason,file,line); break; - case 24: - throw notImplemented(reason,file,line); break; - case 25: - throw cacheLimitExceeded(reason,file,line); break; - case 26: - throw tableLimitExceeded(reason,file,line); break; - case 27: - throw iterationNotProgressing(reason,file,line); break; - case 28: - throw jacobiansNotImprovingSolution(reason,file,line); break; - case 29: - cannotReachToleranceInF(reason,file,line); break; - case 30: - throw cannotReachToleranceInX(reason,file,line); break; - case 31: - throw cannotReachToleranceInGradient(reason,file,line); break; - case 32: - throw endOfFile(reason,file,line); break; - default: //Corresponds to GSL_ERRNO=-1. - throw error(reason,file,line); - } + } - } - - } //ends void handler(const char*, const char*, int, int); +} //ends void handler(const char*, const char*, int, int); }
--- a/src/func.cpp +++ b/src/func.cpp @@ -6,107 +6,107 @@ namespace kwantix{ - //The static variables... - double realfunc::eps = 0; - double realfunc::sqrteps = 0; - double realfunc::root3eps = 0; - bool realfunc::initialised = false; +//The static variables... +double realfunc::eps = 0; +double realfunc::sqrteps = 0; +double realfunc::root3eps = 0; +bool realfunc::initialised = false; - //Static variables - point gsl_function_wrapper::x; - size_t gsl_function_wrapper::index = 1; - realfunc gsl_function_wrapper::myfunc(0); - gsl_function* gsl_function_wrapper::f = 0; +//Static variables +point gsl_function_wrapper::x; +size_t gsl_function_wrapper::index = 1; +realfunc gsl_function_wrapper::myfunc(0); +gsl_function* gsl_function_wrapper::f = 0; - //******************* Wrapper functions ***************************** - gsl_function_wrapper::gsl_function_wrapper( const realfunc &thefunc, - point p, size_t idx){ - myfunc = thefunc; - x = p; - index = idx; - f -> function = &takemyaddress; - f -> params = 0; +//******************* Wrapper functions ***************************** +gsl_function_wrapper::gsl_function_wrapper( const realfunc &thefunc, + point p, size_t idx){ + myfunc = thefunc; + x = p; + index = idx; + f -> function = &takemyaddress; + f -> params = 0; - } +} - void gsl_function_wrapper::set_params(const realfunc &thefunc, - point p, size_t idx){ - myfunc = thefunc; - x = p; - index = idx; - f -> function = &takemyaddress; - f -> params = 0; - } +void gsl_function_wrapper::set_params(const realfunc &thefunc, + point p, size_t idx){ + myfunc = thefunc; + x = p; + index = idx; + f -> function = &takemyaddress; + f -> params = 0; +} - gsl_function* gsl_function_wrapper::get_gsl_function() const{ - return f; - } +gsl_function* gsl_function_wrapper::get_gsl_function() const{ + return f; +} + +double gsl_function_wrapper::takemyaddress(double xi, void* nothing){ + x(index) = xi; + nothing = 0; + return myfunc(x); +} - double gsl_function_wrapper::takemyaddress(double xi, void* nothing){ - x(index) = xi; - nothing = 0; - return myfunc(x); +// **************** realfunc functions ********************************* +realfunc::realfunc() : myfunc(0){ + if(!initialised){ + eps = std::numeric_limits<double>::epsilon(); + sqrteps = sqrt(eps); + root3eps = pow(eps,1/3.0); + initialised = true; } +} - // **************** realfunc functions ********************************* - realfunc::realfunc() : myfunc(0){ - if(!initialised){ - eps = std::numeric_limits<double>::epsilon(); - sqrteps = sqrt(eps); - root3eps = pow(eps,1/3.0); - initialised = true; - } +realfunc::realfunc( double(*f)(const point&)) : myfunc(f) { + if(!initialised){ + eps = std::numeric_limits<double>::epsilon(); + sqrteps = sqrt(eps); + root3eps = pow(eps,1/3.0); + initialised = true; + } +} + +void realfunc::set_function_ptr(double (f_in)(const point &p)){ + myfunc = f_in; +} + +double realfunc::operator()(const point& p) const{ + return at(p); +} +double realfunc::at(const point& p) const{ + if(myfunc == 0){ + throw no_init(__LINE__,__FILE__); } - realfunc::realfunc( double(*f)(const point&)) : myfunc(f) { - if(!initialised){ - eps = std::numeric_limits<double>::epsilon(); - sqrteps = sqrt(eps); - root3eps = pow(eps,1/3.0); - initialised = true; - } - } + return myfunc(p); +} - void realfunc::set_function_ptr(double (f_in)(const point &p)){ - myfunc = f_in; - } +double realfunc::d(const point& p, size_t k) const{ + gsl_function_wrapper gfw(*this,p,k); + double result, abserror; + double x = p(1); + double typx = (1 > log(x) ? 1 : log(x)); + double h = sqrteps*( fabs(x) > typx ? fabs(x) : typx); - double realfunc::operator()(const point& p) const{ - return at(p); - } - double realfunc::at(const point& p) const{ - if(myfunc == 0){ - throw no_init(__LINE__,__FILE__); - } - - return myfunc(p); - } + gsl_deriv_central(gfw.get_gsl_function(), x, h, &result, &abserror); + return result; +} - double realfunc::d(const point& p, size_t k) const{ - gsl_function_wrapper gfw(*this,p,k); - double result, abserror; - double x = p(1); - double typx = (1 > log(x) ? 1 : log(x)); - double h = sqrteps*( fabs(x) > typx ? fabs(x) : typx); - - gsl_deriv_central(gfw.get_gsl_function(), x, h, &result, &abserror); - return result; - } +double realfunc::d2(const point& p, size_t k1, size_t k2) const{ + //FIXME + //Figure this out later. + k1 = k2; + p.size(); + return 0; +} - double realfunc::d2(const point& p, size_t k1, size_t k2) const{ - //FIXME - //Figure this out later. - k1 = k2; - p.size(); - return 0; - } - - kwantix::badArgument - realfunc::no_init(int line, string file) const{ - kwantix::badArgument exc; - exc.line = line; - exc.file = file; - exc.reason = "Did not assign a function pointer to a realfunc object."; - return exc; - } +kwantix::badArgument +realfunc::no_init(int line, string file) const{ + kwantix::badArgument exc; + exc.line = line; + exc.file = file; + exc.reason = "Did not assign a function pointer to a realfunc object."; + return exc; } +}
--- a/src/include/bvp.hpp +++ b/src/include/bvp.hpp @@ -19,204 +19,204 @@ namespace kwantix{ - using std::map; - using std::set; - using std::string; - using kwantix::point; - using kwantix::vector; - using kwantix::matrix; - using boost::shared_ptr; +using std::map; +using std::set; +using std::string; +using kwantix::point; +using kwantix::vector; +using kwantix::matrix; +using boost::shared_ptr; - /*! \brief A domain is a bunch of interior points, a bunch of boundary points, - * and a normal vector at each boundary point. - */ - class domain{ - public: - ///Allocate empty domain of given dimension. - domain(size_t dimension); - ///Allocate domain given a boundary, an interior, and a set of normals. - domain(size_t dimension, set<point> intr, - set<point> bdry, map<point, vector> ns); +/*! \brief A domain is a bunch of interior points, a bunch of boundary points, + * and a normal vector at each boundary point. + */ +class domain{ +public: + ///Allocate empty domain of given dimension. + domain(size_t dimension); + ///Allocate domain given a boundary, an interior, and a set of normals. + domain(size_t dimension, set<point> intr, + set<point> bdry, map<point, vector> ns); - /*!\brief Construct the domain from filenames. - * - * Each file containing the desired points as a matrix structure - * where the number of columns must be the dimensionality of the - * points and each row is one point. - * - * The normals are given by a matrix with number of columns twice - * that of the interior and boundary, with the first part of the - * columns giving the basepoint and the second part of the columns - * giving the normal itself. - * - * A boundary or interior matrix can be empty, but the normals - * matrix can't be empty if the boundary isn't. - * \param intr - Filename holding the interior of the domain. - * \param bdry - Filename holding the boundary of the domain. - * \param ns - Filename holding the normals of the domain. - */ - domain(string intr, string bdry, string ns); + /*!\brief Construct the domain from filenames. + * + * Each file containing the desired points as a matrix structure + * where the number of columns must be the dimensionality of the + * points and each row is one point. + * + * The normals are given by a matrix with number of columns twice + * that of the interior and boundary, with the first part of the + * columns giving the basepoint and the second part of the columns + * giving the normal itself. + * + * A boundary or interior matrix can be empty, but the normals + * matrix can't be empty if the boundary isn't. + * \param intr - Filename holding the interior of the domain. + * \param bdry - Filename holding the boundary of the domain. + * \param ns - Filename holding the normals of the domain. + */ + domain(string intr, string bdry, string ns); - ///This clears any data already in the domain and sets the dimension. - void set_dimension(size_t dimension); + ///This clears any data already in the domain and sets the dimension. + void set_dimension(size_t dimension); - //Add information to the domain - ///Add a set of points to the interior of the domain. - void add_to_interior(const set<point> &intr); - ///Add a point to the interior of the domain. - void add_to_interior(const point &intr); - ///Add a set of points to the boundary of the domain. - void add_to_boundary(const set<point> &bdry); - ///Add a point to the boundary of the domain. - void add_to_boundary(const point &bdry); - /*! \brief Add a set of normals to the domain. - * - * Every normal added to the domain must be attached to a boundary - * point already in the domain. - * \param ns - A set of normals to add. - */ - void add_to_normals(const map<point, vector> &ns); - /*! \brief Add a normal to the domain - * Every normal added to the domain must be attached to a boundary - * point already in the domain. - * \param bdry - The boundary point where to attach this normal. - * \param n - A normal to add. - */ - void add_to_normals(const point &bdry, const vector &n); - - //Read that information - /// Get the domain's dimension. - size_t get_dimension() const; - /// Get the interior. - const set<point>& get_interior() const; - /// Get the boundary. - const set<point>& get_boundary() const; - /// Get the normals. - const map<point, vector>& get_normals() const; - - ///Is point in this domain, whether interior or boundary? - bool contains(const point& p) const; - - virtual ~domain(); - private: - ///Can't create a domain without at least specifying its - ///dimension. - domain(); - ///Domain's dimension - size_t dim; - /*! \name The domain's ingredients - */ - //@{ - set<point> interior; - set<point> boundary; - map<point, vector> normals; - //@} - }; - - /*! \brief A boundary value problem + //Add information to the domain + ///Add a set of points to the interior of the domain. + void add_to_interior(const set<point> &intr); + ///Add a point to the interior of the domain. + void add_to_interior(const point &intr); + ///Add a set of points to the boundary of the domain. + void add_to_boundary(const set<point> &bdry); + ///Add a point to the boundary of the domain. + void add_to_boundary(const point &bdry); + /*! \brief Add a set of normals to the domain. * - * A boundary value problem is a domain \f$\Omega\f$, a differential - * operator \f$\mathcal{L}\f$ on \f$\Omega\f$, a boundary - * differential operator \f$\mathcal{B}\f$ on \f$\partial\Omega\f$, - * and the right hand side values of the equations. Think - * - \f[ - \begin{cases} - \mathcal{L}u = f &\text{on } \Omega \\ - \mathcal{B}u = g &\text{on } \partial\Omega. - \end{cases} - \f] - */ - class BVP{ - public: - /*! \brief Create a boundary value problem. - * - * Given a domain, boundary and interior operators, and values - * that those operators must take on domain and interior, create a - * BVP. - * \param O - A shared_ptr to the domain. - * \param L_in - A shared_ptr to the interior operator. - * \param B_in - A shared_ptr to the boundary operator. - * \param f_in - An std::map or a realfunc giving the values that the - * interior operator must take. - * \param g_in - An std::map or a realfunc giving the values that the - * boundary operator must take. - */ - template<typename F, typename G> - BVP(shared_ptr<const domain> O, - shared_ptr<const diff_op> L_in, - shared_ptr<const bdry_diff_op> B_in, - const F& f_in, - const G& g_in) - { - Omega = O; - L = L_in; - B = B_in; - set_f(f_in); - set_g(g_in); - } + * Every normal added to the domain must be attached to a boundary + * point already in the domain. + * \param ns - A set of normals to add. + */ + void add_to_normals(const map<point, vector> &ns); + /*! \brief Add a normal to the domain + * Every normal added to the domain must be attached to a boundary + * point already in the domain. + * \param bdry - The boundary point where to attach this normal. + * \param n - A normal to add. + */ + void add_to_normals(const point &bdry, const vector &n); - virtual ~BVP() {}; + //Read that information + /// Get the domain's dimension. + size_t get_dimension() const; + /// Get the interior. + const set<point>& get_interior() const; + /// Get the boundary. + const set<point>& get_boundary() const; + /// Get the normals. + const map<point, vector>& get_normals() const; + + ///Is point in this domain, whether interior or boundary? + bool contains(const point& p) const; + + virtual ~domain(); +private: + ///Can't create a domain without at least specifying its + ///dimension. + domain(); + ///Domain's dimension + size_t dim; + /*! \name The domain's ingredients + */ + //@{ + set<point> interior; + set<point> boundary; + map<point, vector> normals; + //@} +}; - ///Get the domain of this BVP. - shared_ptr<const domain> get_domain() const; - ///Get the interior operator of this BVP. - shared_ptr<const diff_op> get_diff_op() const; - ///Get the boundary operator of this BVP. - shared_ptr<const bdry_diff_op> get_bdry_diff_op() const; - ///Get the interior values of this BVP. - const map<point, double>& get_f() const; - ///Get the boundary values of this BVP. - const map<point, double>& get_g() const; +/*! \brief A boundary value problem + * + * A boundary value problem is a domain \f$\Omega\f$, a differential + * operator \f$\mathcal{L}\f$ on \f$\Omega\f$, a boundary + * differential operator \f$\mathcal{B}\f$ on \f$\partial\Omega\f$, + * and the right hand side values of the equations. Think + * + \f[ + \begin{cases} + \mathcal{L}u = f &\text{on } \Omega \\ + \mathcal{B}u = g &\text{on } \partial\Omega. + \end{cases} + \f] +*/ +class BVP{ +public: + /*! \brief Create a boundary value problem. + * + * Given a domain, boundary and interior operators, and values + * that those operators must take on domain and interior, create a + * BVP. + * \param O - A shared_ptr to the domain. + * \param L_in - A shared_ptr to the interior operator. + * \param B_in - A shared_ptr to the boundary operator. + * \param f_in - An std::map or a realfunc giving the values that the + * interior operator must take. + * \param g_in - An std::map or a realfunc giving the values that the + * boundary operator must take. + */ + template<typename F, typename G> + BVP(shared_ptr<const domain> O, + shared_ptr<const diff_op> L_in, + shared_ptr<const bdry_diff_op> B_in, + const F& f_in, + const G& g_in) + { + Omega = O; + L = L_in; + B = B_in; + set_f(f_in); + set_g(g_in); + } - ///Change the interior values of this BVP. - void set_f(const realfunc &f_in); - ///Change the boundary values of this BVP. - void set_g(const realfunc &g_in); + virtual ~BVP() {}; + + ///Get the domain of this BVP. + shared_ptr<const domain> get_domain() const; + ///Get the interior operator of this BVP. + shared_ptr<const diff_op> get_diff_op() const; + ///Get the boundary operator of this BVP. + shared_ptr<const bdry_diff_op> get_bdry_diff_op() const; + ///Get the interior values of this BVP. + const map<point, double>& get_f() const; + ///Get the boundary values of this BVP. + const map<point, double>& get_g() const; + + ///Change the interior values of this BVP. + void set_f(const realfunc &f_in); + ///Change the boundary values of this BVP. + void set_g(const realfunc &g_in); - /*! \brief Change or set interior values of this BVP. - * - * Given an std::map of interior values, this function sets or - * changes the corresponding interior values in the BVP. - * \param f_in - An std::map assigning scalar values to interior - * points. - */ - void set_f(const map<point, double>& f_in); - /*! \brief Change or set interior values of this BVP. - * - * Given an std::map of boundary values, this function sets or - * changes the corresponding boundary values in the BVP. - * \param g_in - An std::map assigning scalar values to boundary - * points. - */ - void set_g(const map<point, double>& g_in); + /*! \brief Change or set interior values of this BVP. + * + * Given an std::map of interior values, this function sets or + * changes the corresponding interior values in the BVP. + * \param f_in - An std::map assigning scalar values to interior + * points. + */ + void set_f(const map<point, double>& f_in); + /*! \brief Change or set interior values of this BVP. + * + * Given an std::map of boundary values, this function sets or + * changes the corresponding boundary values in the BVP. + * \param g_in - An std::map assigning scalar values to boundary + * points. + */ + void set_g(const map<point, double>& g_in); - private: - ///No copying allowed! - BVP(const BVP&){}; - shared_ptr<const domain> Omega; - shared_ptr<const diff_op> L; - shared_ptr<const bdry_diff_op> B; - map<point, double> f; - map<point, double> g; - }; +private: + ///No copying allowed! + BVP(const BVP&){}; + shared_ptr<const domain> Omega; + shared_ptr<const diff_op> L; + shared_ptr<const bdry_diff_op> B; + map<point, double> f; + map<point, double> g; +}; - ///A linear BVP of order at most 2. - class linear_BVP2 : public BVP{ - public: - /// Identical to base class constructor. - template<typename F, typename G> - linear_BVP2(shared_ptr<const domain> O, - shared_ptr<const linear_diff_op2> L_in, - shared_ptr<const bdry_diff_op> B_in, - const F& f_in, - const G& g_in) - : BVP(O, L_in, B_in, f_in, g_in){} +///A linear BVP of order at most 2. +class linear_BVP2 : public BVP{ +public: + /// Identical to base class constructor. + template<typename F, typename G> + linear_BVP2(shared_ptr<const domain> O, + shared_ptr<const linear_diff_op2> L_in, + shared_ptr<const bdry_diff_op> B_in, + const F& f_in, + const G& g_in) + : BVP(O, L_in, B_in, f_in, g_in){} - ///Give the interior diff_op. - shared_ptr<const linear_diff_op2> get_linear_diff_op2() const; - }; + ///Give the interior diff_op. + shared_ptr<const linear_diff_op2> get_linear_diff_op2() const; +}; } #endif //__BVP_HPP__
--- a/src/include/ddm.hpp +++ b/src/include/ddm.hpp @@ -18,106 +18,131 @@ namespace kwantix { - using std::map; - using std::set; - using boost::shared_ptr; +using std::map; +using std::set; +using boost::shared_ptr; - /*! \brief A bdry_diff_op for domain decomposition methods. - * - * Boundary differential operators that may be specified in peculiar - * ways depending on the ddm we're working on. - * - * On interior boundary points as specified by the domain - * decomposition, the alternate Bprime operator is applied instead - * of the B operator. - */ - class ddm_bdry_diff_op : public bdry_diff_op - { - public: - ddm_bdry_diff_op(shared_ptr<const bdry_diff_op> B_in, - shared_ptr<const bdry_diff_op> Bprime_in, - const set<point>& ibps); - double at(const realfunc &f, const point &p) const; - double at(const realfunc &f, const point &p, const vector &n) const; - private: - shared_ptr<const bdry_diff_op> B; - shared_ptr<const bdry_diff_op> Bprime; - set<point> intr_bdry_pts; - }; +/*! \brief A bdry_diff_op for domain decomposition methods. + * + * Boundary differential operators that may be specified in peculiar + * ways depending on the ddm we're working on. + * + * On interior boundary points as specified by the domain + * decomposition, the alternate Bprime operator is applied instead + * of the B operator. + */ +class ddm_bdry_diff_op : public bdry_diff_op +{ +public: + ddm_bdry_diff_op(shared_ptr<const bdry_diff_op> B_in, + shared_ptr<const bdry_diff_op> Bprime_in, + const set<point>& ibps); + double at(const realfunc &f, const point &p) const; + double at(const realfunc &f, const point &p, const vector &n) const; +private: + shared_ptr<const bdry_diff_op> B; + shared_ptr<const bdry_diff_op> Bprime; + set<point> intr_bdry_pts; +}; - //************************ ddm ******************************************** +//************************ ddm ******************************************** - //FIXME: derive this from realfunc - ///A generic domain decomposition method - class ddm{ - public: - ///What are the domains of this ddm, and what's the bvp? - ddm(const set<shared_ptr<const domain> >& ds, - shared_ptr<const BVP> thebvp); +//FIXME: derive this from realfunc +///A generic domain decomposition method +class ddm{ +public: + ///What are the domains of this ddm, and what's the bvp? + ddm(const set<shared_ptr<const domain> >& ds, + shared_ptr<const BVP> thebvp); - ///Relative tolerance. Defaults to 1e-5. - void set_tolerance(double tol); - virtual double at(const point& p) const = 0; - double operator()(const point& p) const; + ///Relative tolerance. Defaults to 1e-5. + void set_tolerance(double tol); + virtual double at(const point& p) const = 0; + double operator()(const point& p) const; - virtual ~ddm(); + virtual ~ddm(); - protected: - /// The BVP for this ddm - shared_ptr<const BVP> bvp; - /// The domain decomposition - set<shared_ptr<const domain> > domains; - /// The tolerance until convergence - double tolerance; +protected: + /// The BVP for this ddm + shared_ptr<const BVP> bvp; + /// The domain decomposition + set<shared_ptr<const domain> > domains; + /// The tolerance until convergence + double tolerance; + + ///Each specific ddm has its own way of doing things. + virtual void solve() = 0; +}; + +class overlapping_domain; + +//***************** additive_schwarz_ddm *********************************** - ///Each specific ddm has its own way of doing things. - virtual void solve() = 0; - }; +/*! \brief Additive Schwarz domain decomposition method. + * + * This method described in Li and Chen [2003]. + */ +template <typename RBF> +class additive_schwarz_ddm : public ddm{ +public: + additive_schwarz_ddm(const set<shared_ptr<const domain> >&, + shared_ptr<const linear_BVP2> bvp); + double at(const point& p) const; +private: + /// Each domain has an associated interpolator.n + map<shared_ptr<const overlapping_domain>, + shared_ptr<interpolator<RBF> > > phis; - class overlapping_domain; - - //***************** additive_schwarz_ddm *********************************** - - /*! \brief Additive Schwarz domain decomposition method. + /*! \brief Say which interpolators are at work at given point. * - * This method described in Li and Chen [2003]. + * \param p - The point + * \return A set of pointers to the interpolators governing this point. */ - template <typename RBF> - class additive_schwarz_ddm : public ddm{ - public: - additive_schwarz_ddm(const set<shared_ptr<const domain> >&, - shared_ptr<const linear_BVP2> bvp); - double at(const point& p) const; - private: - /// Each domain has an associated interpolator.n - map<shared_ptr<const overlapping_domain>, - shared_ptr<interpolator<RBF> > > phis; + set<shared_ptr<const interpolator<RBF> > > + which_interps(const point& p) const; - /*! \brief Say which interpolators are at work at given point. - * - * \param p - The point - * \return A set of pointers to the interpolators governing this point. - */ - set<shared_ptr<const interpolator<RBF> > > - which_interps(const point& p) const; + /// Actual algorithm implemented here, called by constructor. + void solve(); - /// Actual algorithm implemented here, called by constructor. - void solve(); - - ///Evaluate using the average of relevant interpolators. - double avg_interp(set<shared_ptr<const interpolator<RBF> > > - relevant_interpolators, - const point& p) const; + ///Evaluate using the average of relevant interpolators. + double avg_interp(set<shared_ptr<const interpolator<RBF> > > + relevant_interpolators, + const point& p) const; - /*! \brief Return a vector of the current estimated evaluated at all - * artificial boundary points. - */ - vector at_all_points() const; + /*! \brief Return a vector of the current estimated evaluated at all + * artificial boundary points. + */ + vector at_all_points() const; - }; +}; - /*! \brief Overlapping domains for the Scwharz DDMs. +/*! \brief Overlapping domains for the Scwharz DDMs. + * + * An overlapping domain's constructor, besides needing the same + * arguments as required for an ordinary domain, also needs to know + * which domains it will overlap with, and on which domain its + * boundary points lie. These two extra pieces of information are + * given by + * + * set<const shared_ptr<overlapping_domain> > + * map<point, shared_ptr<const overlapping_domain> > + * + * respectively. By default, if these arguments are omitted, it is + * assumed that the overlapping domain doesn't overlap with anything + * at all. If a point in the domain's boundary doesn't overlap with + * another domain, then it shouldn't be specified in the map<point, + * overlapping_domain> data structure (or it can quite harmlessly be + * assigned to itself). + * + * Alternatively, the set_overlapper_info(...) function can handle + * all of the details. + */ +class overlapping_domain : public domain{ +public: + /// A domain that in fact doesn't overlap with anything at all. + overlapping_domain(string intr, string bdry, string ns); + /*! \brief Specify with what other domains it overlaps. * * An overlapping domain's constructor, besides needing the same * arguments as required for an ordinary domain, also needs to know @@ -128,74 +153,49 @@ * set<const shared_ptr<overlapping_domain> > * map<point, shared_ptr<const overlapping_domain> > * - * respectively. By default, if these arguments are omitted, it is - * assumed that the overlapping domain doesn't overlap with anything - * at all. If a point in the domain's boundary doesn't overlap with - * another domain, then it shouldn't be specified in the map<point, - * overlapping_domain> data structure (or it can quite harmlessly be - * assigned to itself). - * - * Alternatively, the set_overlapper_info(...) function can handle - * all of the details. + * respectively. */ - class overlapping_domain : public domain{ - public: - /// A domain that in fact doesn't overlap with anything at all. - overlapping_domain(string intr, string bdry, string ns); - /*! \brief Specify with what other domains it overlaps. - * - * An overlapping domain's constructor, besides needing the same - * arguments as required for an ordinary domain, also needs to know - * which domains it will overlap with, and on which domain its - * boundary points lie. These two extra pieces of information are - * given by - * - * set<const shared_ptr<overlapping_domain> > - * map<point, shared_ptr<const overlapping_domain> > - * - * respectively. - */ - overlapping_domain(string intr, string bdry, string ns, - const set<shared_ptr<const overlapping_domain> >& ols, - const map<point, - shared_ptr<const overlapping_domain> >& bdry_asst); + overlapping_domain(string intr, string bdry, string ns, + const set<shared_ptr<const overlapping_domain> >& ols, + const map<point, + shared_ptr<const overlapping_domain> >& bdry_asst); + + ///Allocate an empty non-overlapping domain of a given dimension + overlapping_domain(size_t dimension); + ///Same, but specifying interior, boundary, and normals. + overlapping_domain(size_t dimension, set<point> intr, + set<point> bdry, map<point, vector> ns); + + /// Give the domains with which this domain overlaps. + set<shared_ptr<const overlapping_domain> > get_domains() const; + + /// With what domain this point overlap? If none, return 0. + shared_ptr<const overlapping_domain> which_domain(const point& p) const; + + friend void set_overlapper_info(set<shared_ptr<overlapping_domain> > + domains); - ///Allocate an empty non-overlapping domain of a given dimension - overlapping_domain(size_t dimension); - ///Same, but specifying interior, boundary, and normals. - overlapping_domain(size_t dimension, set<point> intr, - set<point> bdry, map<point, vector> ns); - - /// Give the domains with which this domain overlaps. - set<shared_ptr<const overlapping_domain> > get_domains() const; - - /// With what domain this point overlap? If none, return 0. - shared_ptr<const overlapping_domain> which_domain(const point& p) const; - - friend void set_overlapper_info(set<shared_ptr<overlapping_domain> > - domains); + ///Set overlapper info manually. + void set_overlapper_info(const point& p, + const shared_ptr<overlapping_domain> o); - ///Set overlapper info manually. - void set_overlapper_info(const point& p, - const shared_ptr<overlapping_domain> o); +private: + /// The domains with which this one overlaps. + set<shared_ptr<const overlapping_domain> > overlappers; + //FIXME: make this a multimap + /// The domains with which this point overlaps. + map<point, shared_ptr<const overlapping_domain> > boundary_assignments; +}; - private: - /// The domains with which this one overlaps. - set<shared_ptr<const overlapping_domain> > overlappers; - //FIXME: make this a multimap - /// The domains with which this point overlaps. - map<point, shared_ptr<const overlapping_domain> > boundary_assignments; - }; - - /*! \brief Automate the initialisation of setting domain - * information. - * - * Once domains are defined via ordinary domain constructors, this - * function will grab a set of domains and say which domain overlaps - * with which one. - */ - void set_overlapper_info(set<shared_ptr<overlapping_domain> > - domains); +/*! \brief Automate the initialisation of setting domain + * information. + * + * Once domains are defined via ordinary domain constructors, this + * function will grab a set of domains and say which domain overlaps + * with which one. + */ +void set_overlapper_info(set<shared_ptr<overlapping_domain> > + domains); } #endif // __DDM_HPP__
--- a/src/include/diff_op.hpp +++ b/src/include/diff_op.hpp @@ -14,139 +14,139 @@ #include "func.hpp" namespace kwantix{ - /* - *Class graph: (other leaves may also be present) - * - * diff_op---------- Dirichlet - * / \ \ / - *linear_diff_op diff_op2 bdry_diff_op--- - * \ / \ - * linear_diff_op2 Neumann - * | - * | - * | - * | - * +-- del1 - * | - * +-- del2 - * | - * +-- Laplacian - * | - * ... etc - * - */ +/* + *Class graph: (other leaves may also be present) + * + * diff_op---------- Dirichlet + * / \ \ / + *linear_diff_op diff_op2 bdry_diff_op--- + * \ / \ + * linear_diff_op2 Neumann + * | + * | + * | + * | + * +-- del1 + * | + * +-- del2 + * | + * +-- Laplacian + * | + * ... etc + * + */ - /*! \brief A general differential operator. - * - * All this class is evaluate the operator applied at a specified - * point for a real-valued function that take a vector - * argument. This is a pure abstract class. +/*! \brief A general differential operator. + * + * All this class is evaluate the operator applied at a specified + * point for a real-valued function that take a vector + * argument. This is a pure abstract class. + */ +class diff_op{ +public: + ///Evaluate the function \f$f\f$ at the point \f$p\f$. + /** All this does is call the pure virtual at(realfunc, point) + * function below. */ - class diff_op{ - public: - ///Evaluate the function \f$f\f$ at the point \f$p\f$. - /** All this does is call the pure virtual at(realfunc, point) - * function below. - */ - double operator()(const realfunc &f, const point &p) const; - /** \brief Pure virtual function. - * - * Actual implementation of the differential operator goes here. - */ - virtual double at(const realfunc &f, const point &p) const = 0; - ///Nothing to destroy. - virtual ~diff_op(){}; - }; + double operator()(const realfunc &f, const point &p) const; + /** \brief Pure virtual function. + * + * Actual implementation of the differential operator goes here. + */ + virtual double at(const realfunc &f, const point &p) const = 0; + ///Nothing to destroy. + virtual ~diff_op(){}; +}; - /// A linear diff_op. Also pure abstract class for now. - class linear_diff_op : virtual public diff_op{ }; +/// A linear diff_op. Also pure abstract class for now. +class linear_diff_op : virtual public diff_op{ }; - /// An at most 2nd-order differential operator. - class diff_op2 : virtual public diff_op{ }; +/// An at most 2nd-order differential operator. +class diff_op2 : virtual public diff_op{ }; - /*! \brief The heat, wave, and Laplace's equation use this kind of - * operators: they're both linear and at most 2nd order. - */ - class linear_diff_op2 : public linear_diff_op, public diff_op2{ - public: - double operator()(const realfunc &f, const point &p) const; - }; +/*! \brief The heat, wave, and Laplace's equation use this kind of + * operators: they're both linear and at most 2nd order. + */ +class linear_diff_op2 : public linear_diff_op, public diff_op2{ +public: + double operator()(const realfunc &f, const point &p) const; +}; - /// An operator for specifying boundary conditions - class bdry_diff_op : public diff_op{ - public: +/// An operator for specifying boundary conditions +class bdry_diff_op : public diff_op{ +public: - double operator()(const realfunc &f, const point &p) const; - /// Only calls at(realfunc, point, vector) below. - double operator()(const realfunc &f, const point &p, - const vector &n) const; - virtual double at(const realfunc &f, - const point &p) const =0; - /*! \brief A boundary differential operator may need to know about - * the normal vector at the boundary. - * - * A boundary differential operator must provide evaluation at - * the boundary, or explicitly ignore it. Robin and Neumann - * boundary operators both make use of this normal vector. - * \param f - The function to which the operator is applied. - * \param p - The point at which the operator is evaluated - * \param n - The normal vector for this point on the boundary. - */ - virtual double at(const realfunc &f, const point &p, - const vector &n) const =0; - }; + double operator()(const realfunc &f, const point &p) const; + /// Only calls at(realfunc, point, vector) below. + double operator()(const realfunc &f, const point &p, + const vector &n) const; + virtual double at(const realfunc &f, + const point &p) const =0; + /*! \brief A boundary differential operator may need to know about + * the normal vector at the boundary. + * + * A boundary differential operator must provide evaluation at + * the boundary, or explicitly ignore it. Robin and Neumann + * boundary operators both make use of this normal vector. + * \param f - The function to which the operator is applied. + * \param p - The point at which the operator is evaluated + * \param n - The normal vector for this point on the boundary. + */ + virtual double at(const realfunc &f, const point &p, + const vector &n) const =0; +}; - /// Dirichlet boundary conditions - class dirichlet_op : public bdry_diff_op{ - public: - double at(const realfunc &f, const point &p) const; - double at(const realfunc &f, const point &p, const vector &n) const; - }; +/// Dirichlet boundary conditions +class dirichlet_op : public bdry_diff_op{ +public: + double at(const realfunc &f, const point &p) const; + double at(const realfunc &f, const point &p, const vector &n) const; +}; - /// Neumann boundary conditions - class neumann_op : public bdry_diff_op{ - public: - double at(const realfunc &f, const point &p, const vector &n) const; - private: - double at(const realfunc &f, const point &p) const {return f(p);}; - }; +/// Neumann boundary conditions +class neumann_op : public bdry_diff_op{ +public: + double at(const realfunc &f, const point &p, const vector &n) const; +private: + double at(const realfunc &f, const point &p) const {return f(p);}; +}; - /// Identity operator - class Id_op : public linear_diff_op2{ - double at(const realfunc &f, const point &p) const; - }; +/// Identity operator +class Id_op : public linear_diff_op2{ + double at(const realfunc &f, const point &p) const; +}; - /*! \brief Partial wrt to some direction. - * - * By default (i.e with the overloaded at() function), wrt first - * coordinate. Else, last argument gives the direction to evaluate - * the partial. - */ - class del1 : public linear_diff_op2{ - public: - double at(const realfunc &f, const point &p) const; - double at(const realfunc &f, const point &p, size_t i) const; - }; +/*! \brief Partial wrt to some direction. + * + * By default (i.e with the overloaded at() function), wrt first + * coordinate. Else, last argument gives the direction to evaluate + * the partial. + */ +class del1 : public linear_diff_op2{ +public: + double at(const realfunc &f, const point &p) const; + double at(const realfunc &f, const point &p, size_t i) const; +}; - /*! \brief Second partials wrt some direction. - * - * By default (i.e with the overloaded at() function), wrt first - * coordinate, twice. Else, last two arguments give the directions - * to evaluate the partials. - */ - class del2 : public linear_diff_op2{ - public: - double at(const realfunc &f, const point &p) const; - double at(const realfunc &f, const point &p, size_t i, size_t j) const; - }; +/*! \brief Second partials wrt some direction. + * + * By default (i.e with the overloaded at() function), wrt first + * coordinate, twice. Else, last two arguments give the directions + * to evaluate the partials. + */ +class del2 : public linear_diff_op2{ +public: + double at(const realfunc &f, const point &p) const; + double at(const realfunc &f, const point &p, size_t i, size_t j) const; +}; - ///The Laplacian - class Laplacian : public linear_diff_op2{ - public: - double at(const realfunc &f, const point &p) const; - }; +///The Laplacian +class Laplacian : public linear_diff_op2{ +public: + double at(const realfunc &f, const point &p) const; +}; }
--- a/src/include/error.hpp +++ b/src/include/error.hpp @@ -19,311 +19,311 @@ #include <string> namespace kwantix{ - //Structs to be thrown as exceptions - using std::string; +//Structs to be thrown as exceptions +using std::string; - /// generic error and base struct. /*GSL_FAILURE = -1,*/ - class error{ - public: - std::string reason, file; - int line; - error(){}; - error(string r, string f, int l) : - reason(r), file(f), line(l) {}; +/// generic error and base struct. /*GSL_FAILURE = -1,*/ +class error{ +public: + std::string reason, file; + int line; + error(){}; + error(string r, string f, int l) : + reason(r), file(f), line(l) {}; - }; +}; - struct noConvergence ; - struct badDomain ; - struct badRange ; - struct badPointer ; - struct badArgument ; - struct failure ; - struct failedFactorisation ; - struct failedSanity ; - struct outOfMemory ; - struct badFunction ; - struct runAway ; - struct maxIterations ; - struct divideByZero ; - struct badTolerance ; - struct aboveTolerance ; - struct underflow ; - struct overflow ; - struct lossOfAccuracy ; - struct roundOffError ; - struct inconformantSizes ; - struct matrixNotSquare ; - struct singularityFound ; - struct integralOrSeriesDivergent ; - struct badHardware ; - struct notImplemented ; - struct cacheLimitExceeded ; - struct tableLimitExceeded ; - struct iterationNotProgressing ; - struct jacobiansNotImprovingSolution ; +struct noConvergence ; +struct badDomain ; +struct badRange ; +struct badPointer ; +struct badArgument ; +struct failure ; +struct failedFactorisation ; +struct failedSanity ; +struct outOfMemory ; +struct badFunction ; +struct runAway ; +struct maxIterations ; +struct divideByZero ; +struct badTolerance ; +struct aboveTolerance ; +struct underflow ; +struct overflow ; +struct lossOfAccuracy ; +struct roundOffError ; +struct inconformantSizes ; +struct matrixNotSquare ; +struct singularityFound ; +struct integralOrSeriesDivergent ; +struct badHardware ; +struct notImplemented ; +struct cacheLimitExceeded ; +struct tableLimitExceeded ; +struct iterationNotProgressing ; +struct jacobiansNotImprovingSolution ; - struct cannotReachToleranceInF ; - struct cannotReachToleranceInX ; - struct cannotReachToleranceInGradient ; - struct endOfFile ; +struct cannotReachToleranceInF ; +struct cannotReachToleranceInX ; +struct cannotReachToleranceInGradient ; +struct endOfFile ; - /*! \brief Custom error handler to be used for GSL. - * - * Throw exceptions instead of calling abort(). - * - * Remember to put `gsl_set_error_handler(&errorHandler);' in the - * main() loops when including this header file; otherwise it's - * useless! - */ - void errorHandler(const char * reason, const char * file, - int line, int gsl_errno); +/*! \brief Custom error handler to be used for GSL. + * + * Throw exceptions instead of calling abort(). + * + * Remember to put `gsl_set_error_handler(&errorHandler);' in the + * main() loops when including this header file; otherwise it's + * useless! + */ +void errorHandler(const char * reason, const char * file, + int line, int gsl_errno); } //A few more details about the structs we're throwing as exceptions. namespace kwantix{ - using std::string; +using std::string; - ///GSL_CONTINUE = -2, /* iteration has not converged */ - struct noConvergence : public error { - noConvergence() {}; - noConvergence(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_CONTINUE = -2, /* iteration has not converged */ +struct noConvergence : public error { + noConvergence() {}; + noConvergence(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EDOM = 1, /* input domain error, e.g sqrt(-1) */ - struct badDomain : public error { - badDomain() {}; - badDomain(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EDOM = 1, /* input domain error, e.g sqrt(-1) */ +struct badDomain : public error { + badDomain() {}; + badDomain(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ERANGE = 2, /* output range error, e.g. exp(1e100) */ - struct badRange : public error { - badRange() {}; - badRange(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ERANGE = 2, /* output range error, e.g. exp(1e100) */ +struct badRange : public error { + badRange() {}; + badRange(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EFAULT = 3, /* invalid pointer */ - struct badPointer : public error { - badPointer() {}; - badPointer(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EFAULT = 3, /* invalid pointer */ +struct badPointer : public error { + badPointer() {}; + badPointer(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EINVAL = 4, /* invalid argument supplied by user */ - struct badArgument : public error { - badArgument() {}; - badArgument(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EINVAL = 4, /* invalid argument supplied by user */ +struct badArgument : public error { + badArgument() {}; + badArgument(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EFAILED = 5, /* generic failure */ - struct failure : public error { - failure() {}; - failure(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EFAILED = 5, /* generic failure */ +struct failure : public error { + failure() {}; + failure(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EFACTOR = 6, /* factorization failed */ - struct failedFactorisation : public error { - failedFactorisation() {}; - failedFactorisation(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EFACTOR = 6, /* factorization failed */ +struct failedFactorisation : public error { + failedFactorisation() {}; + failedFactorisation(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ESANITY = 7, /* sanity check failed - shouldn't happen */ - struct failedSanity : public error { - failedSanity() {}; - failedSanity(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ESANITY = 7, /* sanity check failed - shouldn't happen */ +struct failedSanity : public error { + failedSanity() {}; + failedSanity(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ENOMEM = 8, /* malloc failed */ - struct outOfMemory : public error { - outOfMemory() {}; - outOfMemory(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ENOMEM = 8, /* malloc failed */ +struct outOfMemory : public error { + outOfMemory() {}; + outOfMemory(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EBADFUNC = 9, /* problem with user-supplied function */ - struct badFunction : public error { - badFunction() {}; - badFunction(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EBADFUNC = 9, /* problem with user-supplied function */ +struct badFunction : public error { + badFunction() {}; + badFunction(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ERUNAWAY = 10, /* iterative process is out of control */ - struct runAway : public error { - runAway() {}; - runAway(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ERUNAWAY = 10, /* iterative process is out of control */ +struct runAway : public error { + runAway() {}; + runAway(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EMAXITER = 11, /* exceeded max number of iterations */ - struct maxIterations : public error { - maxIterations() {}; - maxIterations(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EMAXITER = 11, /* exceeded max number of iterations */ +struct maxIterations : public error { + maxIterations() {}; + maxIterations(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EZERODIV = 12, /* tried to divide by zero */ - struct divideByZero : public error { - divideByZero() {}; - divideByZero(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EZERODIV = 12, /* tried to divide by zero */ +struct divideByZero : public error { + divideByZero() {}; + divideByZero(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EBADTOL = 13, /* user specified an invalid tolerance */ - struct badTolerance : public error { - badTolerance() {}; - badTolerance(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EBADTOL = 13, /* user specified an invalid tolerance */ +struct badTolerance : public error { + badTolerance() {}; + badTolerance(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ETOL = 14, /* failed to reach the specified tolerance */ - struct aboveTolerance : public error { - aboveTolerance() {}; - aboveTolerance(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ETOL = 14, /* failed to reach the specified tolerance */ +struct aboveTolerance : public error { + aboveTolerance() {}; + aboveTolerance(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EUNDRFLW = 15, /* underflow */ - struct underflow : public error { - underflow() {}; - underflow(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EUNDRFLW = 15, /* underflow */ +struct underflow : public error { + underflow() {}; + underflow(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EOVRFLW = 16, /* overflow */ - struct overflow : public error { - overflow() {}; - overflow(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EOVRFLW = 16, /* overflow */ +struct overflow : public error { + overflow() {}; + overflow(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ELOSS = 17, /* loss of accuracy */ - struct lossOfAccuracy : public error { - lossOfAccuracy() {}; - lossOfAccuracy(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ELOSS = 17, /* loss of accuracy */ +struct lossOfAccuracy : public error { + lossOfAccuracy() {}; + lossOfAccuracy(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EROUND = 18, /* failed because of roundoff error */ - struct roundOffError : public error { - roundOffError() {}; - roundOffError(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EROUND = 18, /* failed because of roundoff error */ +struct roundOffError : public error { + roundOffError() {}; + roundOffError(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EBADLEN = 19, /* matrix, vector lengths are not conformant */ - struct inconformantSizes : public error { - inconformantSizes() {}; - inconformantSizes(string r, string f, int l) : - error(r,f,l) {}; - size_t n_A, m_A, n_B, m_B; - }; +///GSL_EBADLEN = 19, /* matrix, vector lengths are not conformant */ +struct inconformantSizes : public error { + inconformantSizes() {}; + inconformantSizes(string r, string f, int l) : + error(r,f,l) {}; + size_t n_A, m_A, n_B, m_B; +}; - ///GSL_ENOTSQR = 20, /* matrix not square */ - struct matrixNotSquare : public error { - matrixNotSquare() {}; - matrixNotSquare(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ENOTSQR = 20, /* matrix not square */ +struct matrixNotSquare : public error { + matrixNotSquare() {}; + matrixNotSquare(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ESING = 21, /* apparent singularity detected */ - struct singularityFound : public error { - singularityFound() {}; - singularityFound(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ESING = 21, /* apparent singularity detected */ +struct singularityFound : public error { + singularityFound() {}; + singularityFound(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EDIVERGE = 22, /* integral or series is divergent */ - struct integralOrSeriesDivergent : public error { - integralOrSeriesDivergent() {}; - integralOrSeriesDivergent(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EDIVERGE = 22, /* integral or series is divergent */ +struct integralOrSeriesDivergent : public error { + integralOrSeriesDivergent() {}; + integralOrSeriesDivergent(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EUNSUP = 23, /* requested feature is not supported by the hardware */ - struct badHardware : public error { - badHardware() {}; - badHardware(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EUNSUP = 23, /* requested feature is not supported by the hardware */ +struct badHardware : public error { + badHardware() {}; + badHardware(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EUNIMPL = 24, /* requested feature not (yet) implemented */ - struct notImplemented : public error { - notImplemented() {}; - notImplemented(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EUNIMPL = 24, /* requested feature not (yet) implemented */ +struct notImplemented : public error { + notImplemented() {}; + notImplemented(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ECACHE = 25, /* cache limit exceeded */ - struct cacheLimitExceeded : public error { - cacheLimitExceeded() {}; - cacheLimitExceeded(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ECACHE = 25, /* cache limit exceeded */ +struct cacheLimitExceeded : public error { + cacheLimitExceeded() {}; + cacheLimitExceeded(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ETABLE = 26, /* table limit exceeded */ - struct tableLimitExceeded : public error { - tableLimitExceeded() {}; - tableLimitExceeded(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ETABLE = 26, /* table limit exceeded */ +struct tableLimitExceeded : public error { + tableLimitExceeded() {}; + tableLimitExceeded(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ENOPROG = 27, /* iteration is not making progress towards solution */ - struct iterationNotProgressing : public error { - iterationNotProgressing() {}; - iterationNotProgressing(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ENOPROG = 27, /* iteration is not making progress towards solution */ +struct iterationNotProgressing : public error { + iterationNotProgressing() {}; + iterationNotProgressing(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ENOPROGJ = 28, /* jacobian evaluations are not improving solution */ - struct jacobiansNotImprovingSolution : public error { - jacobiansNotImprovingSolution() {}; - jacobiansNotImprovingSolution(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ENOPROGJ = 28, /* jacobian evaluations are not improving solution */ +struct jacobiansNotImprovingSolution : public error { + jacobiansNotImprovingSolution() {}; + jacobiansNotImprovingSolution(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ETOLF = 29, /* cannot reach the specified tolerance in F */ - struct cannotReachToleranceInF : public error { - cannotReachToleranceInF() {}; - cannotReachToleranceInF(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ETOLF = 29, /* cannot reach the specified tolerance in F */ +struct cannotReachToleranceInF : public error { + cannotReachToleranceInF() {}; + cannotReachToleranceInF(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ETOLX = 30, /* cannot reach the specified tolerance in X */ - struct cannotReachToleranceInX : public error { - cannotReachToleranceInX() {}; - cannotReachToleranceInX(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ETOLX = 30, /* cannot reach the specified tolerance in X */ +struct cannotReachToleranceInX : public error { + cannotReachToleranceInX() {}; + cannotReachToleranceInX(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_ETOLG = 31, /* cannot reach the specified tolerance in gradient */ - struct cannotReachToleranceInGradient : public error { - cannotReachToleranceInGradient() {}; - cannotReachToleranceInGradient(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_ETOLG = 31, /* cannot reach the specified tolerance in gradient */ +struct cannotReachToleranceInGradient : public error { + cannotReachToleranceInGradient() {}; + cannotReachToleranceInGradient(string r, string f, int l) : + error(r,f,l) {}; +}; - ///GSL_EOF = 32 /* end of file */ - struct endOfFile : public error { - endOfFile() {}; - endOfFile(string r, string f, int l) : - error(r,f,l) {}; - }; +///GSL_EOF = 32 /* end of file */ +struct endOfFile : public error { + endOfFile() {}; + endOfFile(string r, string f, int l) : + error(r,f,l) {}; +}; - ///Exception for indices out of range. - struct indexOutOfRange : public badArgument{ - size_t i,j,m,n; - indexOutOfRange() {}; - indexOutOfRange(string r, string f, int l) : - badArgument(r,f,l) {}; - }; +///Exception for indices out of range. +struct indexOutOfRange : public badArgument{ + size_t i,j,m,n; + indexOutOfRange() {}; + indexOutOfRange(string r, string f, int l) : + badArgument(r,f,l) {}; +}; } #endif //__ERROR_HPP__
--- a/src/include/func.hpp +++ b/src/include/func.hpp @@ -19,94 +19,94 @@ namespace kwantix{ - /*! \brief A real-valued function from \f$\mathbb{R}^n\f$ to - * \f$\mathbb{R}\f$ . - * - * That is, a realfunc is a real-valued function that takes in - * kwantix::point and returns doubles. This class is meant to be - * generic and any other real-valued function (e.g. a - * radial_basis_function) should derive from this one. - * - * Note that the base class provides derivatives using - * finite-difference approximations from the GSL. Derived classes - * are intended to provide their own exact derivatives instead of - * using the default finite-difference approximations. - * - */ - class realfunc{ - public: - /// Create a useless realfunc, not assigned to anything. - realfunc(); - /// Create a realfunc from a function pointer. - realfunc( double(*f)(const point&)); - /// Destroys nothing, but declared virtual for derived classes. - virtual ~realfunc(){}; +/*! \brief A real-valued function from \f$\mathbb{R}^n\f$ to + * \f$\mathbb{R}\f$ . + * + * That is, a realfunc is a real-valued function that takes in + * kwantix::point and returns doubles. This class is meant to be + * generic and any other real-valued function (e.g. a + * radial_basis_function) should derive from this one. + * + * Note that the base class provides derivatives using + * finite-difference approximations from the GSL. Derived classes + * are intended to provide their own exact derivatives instead of + * using the default finite-difference approximations. + * + */ +class realfunc{ +public: + /// Create a useless realfunc, not assigned to anything. + realfunc(); + /// Create a realfunc from a function pointer. + realfunc( double(*f)(const point&)); + /// Destroys nothing, but declared virtual for derived classes. + virtual ~realfunc(){}; - /// Give a function pointer to this realfunc. - void set_function_ptr(double (f_in)(const point &p)); + /// Give a function pointer to this realfunc. + void set_function_ptr(double (f_in)(const point &p)); - /// Evaluate the function at point - double operator()(const point& p) const; - /// Evaluate the function at point - virtual double at(const point& p) const; + /// Evaluate the function at point + double operator()(const point& p) const; + /// Evaluate the function at point + virtual double at(const point& p) const; - /// First derivative at point \f$x\f$ and in the \f$k\f$th direction. - virtual double d(const point& x, size_t k) const; - /// Second derivative (FIXME: actually implement this) - virtual double d2(const point& x, size_t k1, size_t k2) const ; - protected: - /// Machine epsilon - static double eps; - /// Square root of epsilon - static double sqrteps; - /// Cube root of epsilon - static double root3eps; - /// Have already initialised epsilon above? - static bool initialised ; - private: - /// Pointer to a function that takes kwantix::point and returns double - double (*myfunc)(const point &p); - /// Exception builder - kwantix::badArgument - no_init(int line, string file) const; - }; + /// First derivative at point \f$x\f$ and in the \f$k\f$th direction. + virtual double d(const point& x, size_t k) const; + /// Second derivative (FIXME: actually implement this) + virtual double d2(const point& x, size_t k1, size_t k2) const ; +protected: + /// Machine epsilon + static double eps; + /// Square root of epsilon + static double sqrteps; + /// Cube root of epsilon + static double root3eps; + /// Have already initialised epsilon above? + static bool initialised ; +private: + /// Pointer to a function that takes kwantix::point and returns double + double (*myfunc)(const point &p); + /// Exception builder + kwantix::badArgument + no_init(int line, string file) const; +}; - ///A function wrapper for calling GSL derivatives. - class gsl_function_wrapper{ - public: - /*! \brief Wraps a realfunc for the GSL - * - * Turns a realfunc into a one-dimensional GSL function from - * \f$\mathbb{R}\f$ to \f$\mathbb{R}\f$ by localising the function - * at a point along a given direction. - * - * \param f - The realfunc to wrap. - * \param p - The point on which to localise. - * \param idx - The direction along with to localise. - */ - gsl_function_wrapper(const realfunc &f, point p, size_t idx); - /// Same as the above constructor - void set_params( const realfunc &f, point p, size_t idx); - /// Pointer to the GSL-style function - gsl_function* get_gsl_function() const; +///A function wrapper for calling GSL derivatives. +class gsl_function_wrapper{ +public: + /*! \brief Wraps a realfunc for the GSL + * + * Turns a realfunc into a one-dimensional GSL function from + * \f$\mathbb{R}\f$ to \f$\mathbb{R}\f$ by localising the function + * at a point along a given direction. + * + * \param f - The realfunc to wrap. + * \param p - The point on which to localise. + * \param idx - The direction along with to localise. + */ + gsl_function_wrapper(const realfunc &f, point p, size_t idx); + /// Same as the above constructor + void set_params( const realfunc &f, point p, size_t idx); + /// Pointer to the GSL-style function + gsl_function* get_gsl_function() const; - /// The actual function that's passed to the GSL - static double takemyaddress(double xi, void* nothing); - private: - /// No nontrivial construction! - gsl_function_wrapper(); + /// The actual function that's passed to the GSL + static double takemyaddress(double xi, void* nothing); +private: + /// No nontrivial construction! + gsl_function_wrapper(); - /// The point at which the function is localised - static point x; - /// The direction along which the function is localised - static size_t index; - /// The function being localised - static realfunc myfunc; + /// The point at which the function is localised + static point x; + /// The direction along which the function is localised + static size_t index; + /// The function being localised + static realfunc myfunc; - /// The GSL function formed with the above data. - static gsl_function* f; - }; + /// The GSL function formed with the above data. + static gsl_function* f; +}; } #endif //__FUNC_HPP__
--- a/src/include/interp_values.hpp +++ b/src/include/interp_values.hpp @@ -5,155 +5,155 @@ #include <map> namespace kwantix{ - using std::map; +using std::map; - class bdry_values; +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; +//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<typename RBF> friend class interpolator; - ///Only constructor available for interpolator class - normals(size_t rbfs_hash_in, const map<point, kwantix::vector>& - normals_in, size_t n_in); - ///Does nothing - normals(){}; - ///Exception builder - badArgument not_initted() 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<typename RBF> friend class interpolator; + ///Only constructor available for interpolator class + normals(size_t rbfs_hash_in, const map<point, kwantix::vector>& + 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; + ///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; + ///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 - kwantix::matrix N; - }; + ///Store normals in a matrix + kwantix::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; +/// 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; - //@} + 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); - //@} + /*! \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 kwantix::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) {}; +protected: + /// Interpolator creates interp_values with this + interp_values(const kwantix::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(); + /// No construction! + interp_values(); - protected: - ///Exception builder - badArgument different_rbfs(int line, string file) const; +protected: + ///Exception builder + badArgument different_rbfs(int line, string file) const; - ///Interpolated values go here - kwantix::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; + ///Interpolated values go here + kwantix::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<typename RBF> friend class interpolator; - }; + template<typename RBF> 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 kwantix::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) {}; - }; +/// 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 kwantix::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) {}; +}; - /// For referring to values purely on the interior of an interpolator. - class intr_values : public interp_values{ - public: - ///Empty initialisation - intr_values(){}; - ///Initialise the interior values with interp_values - intr_values(const interp_values& v_in); - private: - ///Same as for parent class - intr_values(const kwantix::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) {}; - }; +/// For referring to values purely on the interior of an interpolator. +class intr_values : public interp_values{ +public: + ///Empty initialisation + intr_values(){}; + ///Initialise the interior values with interp_values + intr_values(const interp_values& v_in); +private: + ///Same as for parent class + intr_values(const kwantix::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); - //@} +/*! \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__
--- a/src/include/interpolator.hpp +++ b/src/include/interpolator.hpp @@ -22,269 +22,269 @@ #include "interp_values.hpp" namespace kwantix{ - using std::map; - using boost::shared_ptr; +using std::map; +using boost::shared_ptr; - //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. - * - * An interpolator \f$ \displaystyle u(x) = \sum_{\xi \in \Xi} - * \lambda_\xi \varphi_\xi(x) \f$, where the \f$\varphi_\xi\f$ are - * RBFs centred at \f$\xi\f$, and \f$\Xi\f$ is a set of points on - * the interior and boundary of some domain \f$\Omega\f$. - * - * The interpolator can be used in various ways, initialised by - * either a BVP::linear_bvp2 or by interpolation data - * (i.e. std::map< kwantix::point, double>). Once an interpolator is - * initialised, it becomes a bona fide kwantix::realfunc, so that its - * evaluations and derivatives become available. For certain - * problems, it is convenient to be able to evaluate the - * interpolator or its derivatives at all points of the domain. It - * is also useful to be able to interpolate again on the same domain - * or with a slightly modified BVP. Interfaces to both of these - * functions are provided. +//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. + * + * An interpolator \f$ \displaystyle u(x) = \sum_{\xi \in \Xi} + * \lambda_\xi \varphi_\xi(x) \f$, where the \f$\varphi_\xi\f$ are + * RBFs centred at \f$\xi\f$, and \f$\Xi\f$ is a set of points on + * the interior and boundary of some domain \f$\Omega\f$. + * + * The interpolator can be used in various ways, initialised by + * either a BVP::linear_bvp2 or by interpolation data + * (i.e. std::map< kwantix::point, double>). Once an interpolator is + * initialised, it becomes a bona fide kwantix::realfunc, so that its + * evaluations and derivatives become available. For certain + * problems, it is convenient to be able to evaluate the + * interpolator or its derivatives at all points of the domain. It + * is also useful to be able to interpolate again on the same domain + * or with a slightly modified BVP. Interfaces to both of these + * functions are provided. + */ +template<typename RBF> +class interpolator : public realfunc{ +public: + + /** @name Constructors + * + * Constructors that take interpolation data perform the + * interpolation as part of the initialisation, which includes + * factoring a matrix */ - template<typename RBF> - class interpolator : public realfunc{ - public: - - /** @name Constructors - * - * Constructors that take interpolation data perform the - * interpolation as part of the initialisation, which includes - * factoring a matrix - */ - //@{ - ///Does not initialise the interpolator. - interpolator(); + //@{ + ///Does not initialise the interpolator. + interpolator(); - ///Interpolate given a BVP - interpolator(shared_ptr<linear_BVP2> bvp); + ///Interpolate given a BVP + interpolator(shared_ptr<linear_BVP2> bvp); - ///Interpolate given some data points and the value at those points - interpolator(const map<point, double>& Xi); + ///Interpolate given some data points and the value at those points + interpolator(const map<point, double>& Xi); - /** \short Interpolate given a domain, data points and the value - * at those points - * - * Must provide domain information. The values of Xi must match - * points on the given domain Omega. The reason why the domain - * might be relevant is in case it is important to specify which - * points are on the interior and which are on the boundary, so - * that it's possible to interpolate again by changing the values - * on either the interior or the boundary. - */ - interpolator(shared_ptr<domain> Omega, const map<point, double>& Xi); - //@} + /** \short Interpolate given a domain, data points and the value + * at those points + * + * Must provide domain information. The values of Xi must match + * points on the given domain Omega. The reason why the domain + * might be relevant is in case it is important to specify which + * points are on the interior and which are on the boundary, so + * that it's possible to interpolate again by changing the values + * on either the interior or the boundary. + */ + interpolator(shared_ptr<domain> Omega, const map<point, double>& Xi); + //@} - /** @name Interpolation - * - * Interpolate again either given new data or a different BVP. - */ - //@{ - void interpolate(const map<point, double>& Xi); - void interpolate(shared_ptr<linear_BVP2> bvp); + /** @name Interpolation + * + * Interpolate again either given new data or a different BVP. + */ + //@{ + void interpolate(const map<point, double>& Xi); + void interpolate(shared_ptr<linear_BVP2> bvp); - /** \short 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); - //@} + /** \short 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 Pointwise evaluation and derivatives - */ - //@{ - ///Point evaluation - double operator()(const point& p) const; - ///Point evaluation - double at(const point& p) const; - /// Pointwise first derivative - double d(const point& p, size_t k) const; - /// Pointwise second derivatives - double d2(const point &p, size_t k1, size_t k2) const; - //@} + /** @name Pointwise evaluation and derivatives + */ + //@{ + ///Point evaluation + double operator()(const point& p) const; + ///Point evaluation + double at(const point& p) const; + /// Pointwise first derivative + double d(const point& p, size_t k) const; + /// Pointwise second derivatives + double d2(const point &p, size_t k1, size_t k2) const; + //@} - /** @name Full domain evaluation and derivatives - */ - //@{ - ///Full domain evaluation - interp_values operator()() const; - ///Full domain evaluation - interp_values at() const; - /// Full domain first derivative - interp_values d(size_t k) const; - /// Full domain second derivatives - interp_values d2(size_t k1, size_t k2) const; - //@} + /** @name Full domain evaluation and derivatives + */ + //@{ + ///Full domain evaluation + interp_values operator()() const; + ///Full domain evaluation + interp_values at() const; + /// Full domain first derivative + interp_values d(size_t k) const; + /// Full domain second derivatives + interp_values d2(size_t k1, size_t k2) const; + //@} - /** @name Precomputation of RBFs - * - * The following functions will precompute upon request all the - * interpolator's RBFs at all points of the domain. This may speed - * up some computations at the expense of increasing memory - * usage. - */ - //@{ - ///Precompute all the function evaluation of the RBFs. - void precompute_ev(); - ///Precompute all relevant first derivatives. - void precompute_d1() ; - ///Precompute all relevant second derivatives. - void precompute_d2(); - //@} + /** @name Precomputation of RBFs + * + * The following functions will precompute upon request all the + * interpolator's RBFs at all points of the domain. This may speed + * up some computations at the expense of increasing memory + * usage. + */ + //@{ + ///Precompute all the function evaluation of the RBFs. + void precompute_ev(); + ///Precompute all relevant first derivatives. + void precompute_d1() ; + ///Precompute all relevant second derivatives. + void precompute_d2(); + //@} - /** @name Partial redefinitions - * - * These functions allow for partial redefinition of the BVP as - * equired for the additive Schwartz domain decomposition method, - * and for other methods. They do not factor a matrix again. - */ - //@{ - void set_f(const intr_values &f); - void set_g(const bdry_values &g); - void set_f(const realfunc &f); - void set_g(const realfunc &g); - void set_f(const map<point, double>& f); - void set_g(const map<point, double>& g); - //@} + /** @name Partial redefinitions + * + * These functions allow for partial redefinition of the BVP as + * equired for the additive Schwartz domain decomposition method, + * and for other methods. They do not factor a matrix again. + */ + //@{ + void set_f(const intr_values &f); + void set_g(const bdry_values &g); + void set_f(const realfunc &f); + void set_g(const realfunc &g); + void set_f(const map<point, double>& f); + 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 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 - * linear operations and work by manipulating interpolator - * coefficients. - */ - //@{ - /// Needs two operators on the same domain. - interpolator<RBF> operator+(const interpolator<RBF>& u) const; - /// Needs two operators on the same domain. - interpolator<RBF> operator-(const interpolator<RBF>& u) const; - /// Returns a scaled interpolator. - interpolator<RBF> operator*(double a) const; - /// Returns a scaled interpolator. - interpolator<RBF> operator/(double a) const; - //@} + /** @name Linear arithmetic operators + * + * These functions return a new interpolator. They are pointwise + * linear operations and work by manipulating interpolator + * coefficients. + */ + //@{ + /// Needs two operators on the same domain. + interpolator<RBF> operator+(const interpolator<RBF>& u) const; + /// Needs two operators on the same domain. + interpolator<RBF> operator-(const interpolator<RBF>& u) const; + /// Returns a scaled interpolator. + interpolator<RBF> operator*(double a) const; + /// Returns a scaled interpolator. + interpolator<RBF> operator/(double a) const; + //@} - /** \brief Provide a unit normals object associated to the - * interpolator. - */ - normals get_normals() const; + /** \brief Provide a unit normals object associated to the + * interpolator. + */ + normals get_normals() const; - /** \brief This will put into an ostream a kwantix::matrix where - * the first \f$n-1\f$ columns will be all the points where the - * interpolator is defined and the \f$n^{th}\f$ column is the - * values that the interpolator takes on that point. - */ - void into_os(std::ostream& os) const; + /** \brief This will put into an ostream a kwantix::matrix where + * the first \f$n-1\f$ columns will be all the points where the + * interpolator is defined and the \f$n^{th}\f$ column is the + * values that the interpolator takes on that point. + */ + void into_os(std::ostream& os) const; - private: - ///Once the matrix is defined, this function inverts it. - void computecoeffs(bool rhs_defined = false); +private: + ///Once the matrix is defined, this function inverts it. + void computecoeffs(bool rhs_defined = false); - ///Perform the actual interpolation. - void init(shared_ptr<linear_BVP2> bvp); + ///Perform the actual interpolation. + void init(shared_ptr<linear_BVP2> bvp); - ///If the coefficients change, need to recompute precomputed values - void recompute_values_vec(); + ///If the coefficients change, need to recompute precomputed values + void recompute_values_vec(); - ///BVP associated to this interpolator - shared_ptr<linear_BVP2> thebvp; + ///BVP associated to this interpolator + shared_ptr<linear_BVP2> thebvp; - ///Number of interior points. - size_t n; - ///Number of boundary points. - size_t m; + ///Number of interior points. + size_t n; + ///Number of boundary points. + size_t m; - ///The matrix to invert. - matrix M; + ///The matrix to invert. + matrix M; - ///Is the interpolator ready for use? - bool initted; + ///Is the interpolator ready for use? + bool initted; - /** @name Exception builders - */ - //@{ - kwantix::badArgument not_initted(int line, string file) const; - kwantix::badArgument not_precomputed(int line, string file) const; - //@} + /** @name Exception builders + */ + //@{ + kwantix::badArgument not_initted(int line, string file) const; + kwantix::badArgument not_precomputed(int line, string file) const; + //@} - ///Coefficients of the RBFs - kwantix::vector coeffs; - ///The RBFS - std::vector<RBF> rbfs; - ///Their hash - size_t rbfs_hash; + ///Coefficients of the RBFs + kwantix::vector coeffs; + ///The RBFS + std::vector<RBF> rbfs; + ///Their hash + size_t rbfs_hash; - ///An RBF hasher - size_t hash_value(const std::vector<RBF>& rbfs_in); + ///An RBF hasher + size_t hash_value(const std::vector<RBF>& rbfs_in); - ///The normals at the domain's boundary - normals nrmls; + ///The normals at the domain's boundary + normals nrmls; - /** \short Precomputed RBFs - * - * For all points on the domain, this stores all the RBFs - * evaluated at those points, as well as the derivatives on where - * the RBFs have been evaluated. This is to speed up successive - * evaluations of the interpolator when the interpolator's domain - * doesn't change. - * - * 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; + /** \short Precomputed RBFs + * + * For all points on the domain, this stores all the RBFs + * evaluated at those points, as well as the derivatives on where + * the RBFs have been evaluated. This is to speed up successive + * evaluations of the interpolator when the interpolator's domain + * doesn't change. + * + * 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; - /** @name Evaluation flags - *Have the RBFs or their derivatives been precomputed? - */ - //@{ - bool ev_precomp; - bool d1_precomp; - bool d2_precomp; - //@} + /** @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>, kwantix::vector > - precomp_values_vec; - /// Convert vector of values to a map - map<point, double> valsvec2map(const kwantix::vector& values) const; - }; + /// 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>, kwantix::vector > + precomp_values_vec; + /// Convert vector of values to a map + map<point, double> valsvec2map(const kwantix::vector& values) const; +}; - /// For comfortable syntax - template <typename RBF> - interpolator<RBF> operator*(double a, const interpolator<RBF>& u) - { - return u*a; - } +/// For comfortable syntax +template <typename RBF> +interpolator<RBF> operator*(double a, const interpolator<RBF>& u) +{ + return u*a; +} - /// Stream insertion - template <typename RBF> - std::ostream& operator<< (std::ostream& os, const interpolator<RBF>& u) - { - u.into_os(os); - return os; - } +/// Stream insertion +template <typename RBF> +std::ostream& operator<< (std::ostream& os, const interpolator<RBF>& u) +{ + u.into_os(os); + return os; +} }
--- a/src/include/linalg.hpp +++ b/src/include/linalg.hpp @@ -22,474 +22,474 @@ namespace kwantix{ - class slice; - class vector; - class vector_view; +class slice; +class vector; +class vector_view; - /// A wrapper class for GNU Scientific Library matrices - class matrix{ - private: - class LUmatrix; +/// A wrapper class for GNU Scientific Library matrices +class matrix{ +private: + class LUmatrix; - public: - //Constructors +public: + //Constructors - /// Allocate 1x1 matrix zero matrix. - matrix(); - /*! \brief Allocate m by n matrix, filled with some value. - \param m - Number of rows. - \param n - Number of columns. - \param fillvalue - Value to fill matrix with, default is 0. - */ - matrix(const size_t m, const size_t n, const double fillvalue = 0) ; - /// Assign a matrix from a GSL matrix pointer. - matrix(gsl_matrix *M); + /// Allocate 1x1 matrix zero matrix. + matrix(); + /*! \brief Allocate m by n matrix, filled with some value. + \param m - Number of rows. + \param n - Number of columns. + \param fillvalue - Value to fill matrix with, default is 0. + */ + matrix(const size_t m, const size_t n, const double fillvalue = 0) ; + /// Assign a matrix from a GSL matrix pointer. + matrix(gsl_matrix *M); - /// Copy constructor - matrix(const matrix& M); - /// Assignment. - matrix operator=(const matrix& M); + /// Copy constructor + matrix(const matrix& M); + /// Assignment. + matrix operator=(const matrix& M); - //Destructors + //Destructors - /// Clears the memory allocated by the GSL - ~matrix(); + /// Clears the memory allocated by the GSL + ~matrix(); - /// Number of decimal digits to output. - size_t precision() const; - /// Set the number of decimal digits to output. - static void set_precision(size_t p); + /// Number of decimal digits to output. + size_t precision() const; + /// Set the number of decimal digits to output. + static void set_precision(size_t p); - /// Number of rows. - size_t rows() const; - /// Number of columns. - size_t cols() const; + /// Number of rows. + size_t rows() const; + /// Number of columns. + size_t cols() const; - /*! \brief Fortran-style parenthetical indexing (hence - * Octave-style too) - * - * Indices start from 1. - * @param i - Row number. - * @param j - Column number. - * @return A reference to the matrix element. - */ - double& operator()(const size_t i, const size_t j); + /*! \brief Fortran-style parenthetical indexing (hence + * Octave-style too) + * + * Indices start from 1. + * @param i - Row number. + * @param j - Column number. + * @return A reference to the matrix element. + */ + double& operator()(const size_t i, const size_t j); - /// Constant version of previous function. - const double& operator()(const size_t i, const size_t j) const; + /// Constant version of previous function. + const double& operator()(const size_t i, const size_t j) const; - /*! \brief Indexing for vectorisation. - * - * The GSL provides limited vectorisation routines which can be - * useful for operating on blocks of matrices all at once. - * @param i - Row to be sliced. - * @param b - Slice range. - * @return A vector_view of the sliced row. - */ - vector_view operator()(const size_t i, const slice &b); - /// Constant version of previous function. - const vector_view operator()(const size_t i, const slice &b) const; + /*! \brief Indexing for vectorisation. + * + * The GSL provides limited vectorisation routines which can be + * useful for operating on blocks of matrices all at once. + * @param i - Row to be sliced. + * @param b - Slice range. + * @return A vector_view of the sliced row. + */ + vector_view operator()(const size_t i, const slice &b); + /// Constant version of previous function. + const vector_view operator()(const size_t i, const slice &b) const; - /*! \brief Indexing for vectorisation. - * - * The GSL provides limited vectorisation routines which can be - * useful for operating on blocks of matrices all at once. - * @param a - Slice range - * @param j - column to be sliced. - * @return A vector_view of the sliced column. - */ - vector_view operator()(const slice &a, const size_t j); - /// Constant version of previous function. - const vector_view operator()(const slice &a, const size_t j) const; + /*! \brief Indexing for vectorisation. + * + * The GSL provides limited vectorisation routines which can be + * useful for operating on blocks of matrices all at once. + * @param a - Slice range + * @param j - column to be sliced. + * @return A vector_view of the sliced column. + */ + vector_view operator()(const slice &a, const size_t j); + /// Constant version of previous function. + const vector_view operator()(const slice &a, const size_t j) const; - //Arithmetic operations - ///Scale the matrix. - matrix operator*(const double a) const; - /// Matrix addition - matrix operator+(const matrix& N) const; - /// Matrix multiplication - matrix operator*(const matrix& N) const; - /// Matrix subtraction - matrix operator-(const matrix& N) const; - /// Matrix vector product (gaxpy operation). - vector operator*(const vector& v) const; //Mv, where v is treated - //as a column vector. + //Arithmetic operations + ///Scale the matrix. + matrix operator*(const double a) const; + /// Matrix addition + matrix operator+(const matrix& N) const; + /// Matrix multiplication + matrix operator*(const matrix& N) const; + /// Matrix subtraction + matrix operator-(const matrix& N) const; + /// Matrix vector product (gaxpy operation). + vector operator*(const vector& v) const; //Mv, where v is treated + //as a column vector. - //More complex operations. + //More complex operations. - /*! \brief Inverse. - * - * This function returns the inverse matrix, computed via an LU - * factorisation. If the matrix is ill-conditioned or the inverse - * does not exist, it will happily return a matrix full of NaN. - * \returns The inverse matrix. - */ - matrix inv() const; - /*! \brief Tranpose. - * - * Returns a copy of the original matrix, transposed. The original - * matrix is not modified. - * \returns Transposed copy of matrix. - */ - matrix T() const; - /*! \brief Trace. - * - * The trace of a square matrix. - * \returns The trace. - */ - double tr() const; - /*! \brief Determinant. - * - * The determinant computed via an LU factorisation - * \returns The determinant. - */ - double det() const; + /*! \brief Inverse. + * + * This function returns the inverse matrix, computed via an LU + * factorisation. If the matrix is ill-conditioned or the inverse + * does not exist, it will happily return a matrix full of NaN. + * \returns The inverse matrix. + */ + matrix inv() const; + /*! \brief Tranpose. + * + * Returns a copy of the original matrix, transposed. The original + * matrix is not modified. + * \returns Transposed copy of matrix. + */ + matrix T() const; + /*! \brief Trace. + * + * The trace of a square matrix. + * \returns The trace. + */ + double tr() const; + /*! \brief Determinant. + * + * The determinant computed via an LU factorisation + * \returns The determinant. + */ + double det() const; - /*! \brief Solves Mv = w for v with LU factorisation. - * - * Given another vector, will return the solution vector v to the - * equation Mv = w. - * \param w - Another vector, must have the same size as number of - * rows in the matrix. - * \returns The solution v to the equation Mv = w. - */ - vector inv(const vector& w) const; + /*! \brief Solves Mv = w for v with LU factorisation. + * + * Given another vector, will return the solution vector v to the + * equation Mv = w. + * \param w - Another vector, must have the same size as number of + * rows in the matrix. + * \returns The solution v to the equation Mv = w. + */ + vector inv(const vector& w) const; - /*! \brief L2 condition number, using SVD. - * - * This function computes the L2 condition number of a matrix, the - * largest singular value divided by the smallest. It directly - * uses the SVD decomposition (as provided by the GSL), and it may - * take a long time for large matrices. - * \returns The L2 condition number. - */ - double cond() const; + /*! \brief L2 condition number, using SVD. + * + * This function computes the L2 condition number of a matrix, the + * largest singular value divided by the smallest. It directly + * uses the SVD decomposition (as provided by the GSL), and it may + * take a long time for large matrices. + * \returns The L2 condition number. + */ + double cond() const; - friend class vector_view; - private: - /*! \brief LU decomposition, pivots in U. - * - * This function will compute the LU factorisation of a matrix, - * where the pivots are in the upper triangular matrix U and L has - * ones along the diagonal (not stored). For efficiency, L and U - * are both internally stored within a single matrix. - * \returns A pointer to the private class LUmatrix. - */ - LUmatrix* LU() const; + friend class vector_view; +private: + /*! \brief LU decomposition, pivots in U. + * + * This function will compute the LU factorisation of a matrix, + * where the pivots are in the upper triangular matrix U and L has + * ones along the diagonal (not stored). For efficiency, L and U + * are both internally stored within a single matrix. + * \returns A pointer to the private class LUmatrix. + */ + LUmatrix* LU() const; - gsl_matrix * A; + gsl_matrix * A; - ///Precision - static size_t precsn; - /// Machine epsilon - static const double eps; + ///Precision + static size_t precsn; + /// Machine epsilon + static const double eps; - // Internal members, all can change, none of them affect the - // actual matrix per se. - mutable LUmatrix* LUptr; - mutable bool LUfactored; + // Internal members, all can change, none of them affect the + // actual matrix per se. + mutable LUmatrix* LUptr; + mutable bool LUfactored; - mutable double condition_number; //L2 condition number, obtained by svd. - mutable gsl_vector* SVD; //Matrix's singular values. - mutable bool SVDfactored; + mutable double condition_number; //L2 condition number, obtained by svd. + mutable gsl_vector* SVD; //Matrix's singular values. + mutable bool SVDfactored; - void SVDfactor() const; + void SVDfactor() const; - /// A private matrix member class for matrices factored in LU form. - class LUmatrix{ - public: - gsl_matrix* matrix_ptr(); - gsl_permutation* perm_ptr(); - int sgn(); - int* sgn_ptr(); + /// A private matrix member class for matrices factored in LU form. + class LUmatrix{ + public: + gsl_matrix* matrix_ptr(); + gsl_permutation* perm_ptr(); + int sgn(); + int* sgn_ptr(); - LUmatrix(gsl_matrix* M); - LUmatrix(const LUmatrix& LU); - LUmatrix operator=(const LUmatrix& LU); - ~LUmatrix(); - private: - // Private interface to the GSL. - gsl_permutation* p; - gsl_matrix* A; - int signum; - LUmatrix(); - };//End LUmatrix - }; //End matrix + LUmatrix(gsl_matrix* M); + LUmatrix(const LUmatrix& LU); + LUmatrix operator=(const LUmatrix& LU); + ~LUmatrix(); + private: + // Private interface to the GSL. + gsl_permutation* p; + gsl_matrix* A; + int signum; + LUmatrix(); + };//End LUmatrix +}; //End matrix - /// A wrapper class for GSL vectors - class vector{ - public: - //Constructors +/// A wrapper class for GSL vectors +class vector{ +public: + //Constructors - ///Allocate zero vector of size one. - vector(); - /*! \brief Allocate GSL vector of size n, filled with fillvalue - * \param n - Size of new vector - * \param fillvalue - Value for filling vector, default is 0. - */ - vector(const size_t n, const double fillvalue = 0); - /// Associate a vector to a GSL vector pointer. - vector(gsl_vector *y); - /// Associate a vector to a constant GSL vector pointer. - vector(const gsl_vector *y); + ///Allocate zero vector of size one. + vector(); + /*! \brief Allocate GSL vector of size n, filled with fillvalue + * \param n - Size of new vector + * \param fillvalue - Value for filling vector, default is 0. + */ + vector(const size_t n, const double fillvalue = 0); + /// Associate a vector to a GSL vector pointer. + vector(gsl_vector *y); + /// Associate a vector to a constant GSL vector pointer. + vector(const gsl_vector *y); - /// Copy contstructor. - vector(const vector& y); - /// Assignment. - vector& operator=(const vector &y); + /// Copy contstructor. + vector(const vector& y); + /// Assignment. + vector& operator=(const vector &y); - //Destructor - /// Clear memory assigned to vector by GSL. - ~vector(); + //Destructor + /// Clear memory assigned to vector by GSL. + ~vector(); - //Number of decimal digits to output. - /// Number of decimal digits to output - size_t precision() const; - /// Set the number of decimal digits to output - static void set_precision(size_t p); + //Number of decimal digits to output. + /// Number of decimal digits to output + size_t precision() const; + /// Set the number of decimal digits to output + static void set_precision(size_t p); - ///Number of elements - size_t size() const; + ///Number of elements + size_t size() const; - /*! \brief Fortran-style parenthetical indexing (hence - * Octave-style too). Indices start at 1. - * \param i - Index number. - * \return A reference to the vector element. - */ - double& operator()(const size_t i); + /*! \brief Fortran-style parenthetical indexing (hence + * Octave-style too). Indices start at 1. + * \param i - Index number. + * \return A reference to the vector element. + */ + double& operator()(const size_t i); - /// Constant version of previous function. - const double& operator()(const size_t i) const; + /// Constant version of previous function. + const double& operator()(const size_t i) const; - /*! \brief Indexing for vectorisation. - * The GSL provides limited vectorisation routines which can be - * useful for operating on blocks of vectors all at once. - * @param v - Slice range. - * @return A vector_view of the slice. - */ - vector_view operator()(const slice &v); - ///Constant version of previous function. - const vector_view operator()(const slice &v) const; + /*! \brief Indexing for vectorisation. + * The GSL provides limited vectorisation routines which can be + * useful for operating on blocks of vectors all at once. + * @param v - Slice range. + * @return A vector_view of the slice. + */ + vector_view operator()(const slice &v); + ///Constant version of previous function. + const vector_view operator()(const slice &v) const; - //Arithmetic operations: - ///Scale the vector. - vector operator*(const double a) const; - ///Scale the vector. - vector operator/(const double a) const; - ///Vector addition. - vector operator+(const vector& w) const; - ///Vector subtration. - vector 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; + //Arithmetic operations: + ///Scale the vector. + vector operator*(const double a) const; + ///Scale the vector. + vector operator/(const double a) const; + ///Vector addition. + vector operator+(const vector& w) const; + ///Vector subtration. + vector 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; - //Comparison + //Comparison - /// Compares vectors elementwise, using machine epsilon. - bool operator==(const vector& w) const; - /*! \brief Lexicographical order, used for putting into STL sets or - * maps. - * - * Lexicographical ordering of vectors. Vectors of smaller - * dimension are smaller. - */ - bool operator<(const vector& w) const; + /// Compares vectors elementwise, using machine epsilon. + bool operator==(const vector& w) const; + /*! \brief Lexicographical order, used for putting into STL sets or + * maps. + * + * Lexicographical ordering of vectors. Vectors of smaller + * dimension are smaller. + */ + bool operator<(const vector& w) const; - //More complex operations - ///Euclidean norm. - double norm() const; + //More complex operations + ///Euclidean norm. + double norm() const; - ///Dot product - double operator%(const vector& w) const; + ///Dot product + double operator%(const vector& w) const; - friend class vector_view; - friend class matrix; - private: - /// Pointer to associated GSL vector. - gsl_vector * x; - /// Output precision. - static size_t precsn; - /// Machine epsilon - static const double eps; - }; + friend class vector_view; + friend class matrix; +private: + /// Pointer to associated GSL vector. + gsl_vector * x; + /// Output precision. + static size_t precsn; + /// Machine epsilon + static const double eps; +}; - /*! \brief A vector that doesn't own its data; rather, points to data owned - * by another vector. - */ - class vector_view : public vector{ - friend class vector; - friend class matrix; - public: - ~vector_view(); - /// Assignment to a vector; will share same data in memory as that - /// other vector. - vector_view& operator=(const vector& w); - /// Assignment, useful for modifying chunks of vectors all at once. - vector_view& operator=(const vector_view& w); - private: - /// vector_views are invisible to the user and can only - /// be created using slices and matrix or vector indexing. - vector_view(); - /// See above. - vector_view(gsl_vector* y, const slice& s); - /// See above. - vector_view(gsl_matrix* A, const slice& a, const size_t j); - /// See above. - vector_view(gsl_matrix* A, const size_t i, const slice& b); - }; +/*! \brief A vector that doesn't own its data; rather, points to data owned + * by another vector. + */ +class vector_view : public vector{ + friend class vector; + friend class matrix; +public: + ~vector_view(); + /// Assignment to a vector; will share same data in memory as that + /// other vector. + vector_view& operator=(const vector& w); + /// Assignment, useful for modifying chunks of vectors all at once. + vector_view& operator=(const vector_view& w); +private: + /// vector_views are invisible to the user and can only + /// be created using slices and matrix or vector indexing. + vector_view(); + /// See above. + vector_view(gsl_vector* y, const slice& s); + /// See above. + vector_view(gsl_matrix* A, const slice& a, const size_t j); + /// See above. + vector_view(gsl_matrix* A, const size_t i, const slice& b); +}; - /*! \brief Vector slices corresponding to GNU Octave ranges. - * - * Vector slices are useful for certain vectorisation routines and - * for referring to rows, columns, or equally spaced chunks of a - * vector or a matrix. Set a begininning and an end for the - * stride, and a step size, and this will be the elements that can - * be indexed by this slice in a vector or matrix. E.g. if the - * beginning is 2, the end 8, and the stride is 2, then this - * refers to the slice [2, 4, 6, 8]. The default stride is 1. - * - * Unlike GNU Octave slices, there is no shorthand for referring - * to the entirety of a column or row of a matrix. Must instead - * create a slice for that purpose as s(1,n) where n is the number - * of rows or columns respectively. +/*! \brief Vector slices corresponding to GNU Octave ranges. + * + * Vector slices are useful for certain vectorisation routines and + * for referring to rows, columns, or equally spaced chunks of a + * vector or a matrix. Set a begininning and an end for the + * stride, and a step size, and this will be the elements that can + * be indexed by this slice in a vector or matrix. E.g. if the + * beginning is 2, the end 8, and the stride is 2, then this + * refers to the slice [2, 4, 6, 8]. The default stride is 1. + * + * Unlike GNU Octave slices, there is no shorthand for referring + * to the entirety of a column or row of a matrix. Must instead + * create a slice for that purpose as s(1,n) where n is the number + * of rows or columns respectively. + */ +class slice{ +public: + /*! \brief Create a vector slice. + * + * \param a - Beginning of the slice. + * \param b - End of the slice. + * \param k - Stride (default is 1). */ - class slice{ - public: - /*! \brief Create a vector slice. - * - * \param a - Beginning of the slice. - * \param b - End of the slice. - * \param k - Stride (default is 1). - */ - slice(size_t a, size_t b, size_t k=1); + slice(size_t a, size_t b, size_t k=1); + + /*! For setting the slice parameters anew. + * + * Indices start from 1. + * \param a - Beginning of the slice. + * \param b - End of the slice. + * \param k - Stride (default is 1). + */ + slice operator()(size_t a, size_t b, size_t k=1); + slice set(size_t a, size_t b, size_t k=1) ; - /*! For setting the slice parameters anew. - * - * Indices start from 1. - * \param a - Beginning of the slice. - * \param b - End of the slice. - * \param k - Stride (default is 1). - */ - slice operator()(size_t a, size_t b, size_t k=1); - slice set(size_t a, size_t b, size_t k=1) ; + size_t begin()const {return beg;}; + size_t end() const{return fin;}; + size_t stride() const{return str;}; - size_t begin()const {return beg;}; - size_t end() const{return fin;}; - size_t stride() const{return str;}; - - private: - /// Empty and private default constructor. - slice(){}; - size_t beg,fin; - size_t str; - }; +private: + /// Empty and private default constructor. + slice(){}; + size_t beg,fin; + size_t str; +}; - /// Useful alias, vectors are also points in space. - typedef vector point; //Useful alias. +/// Useful alias, vectors are also points in space. +typedef vector point; //Useful alias. } namespace kwantix{ //Non-member functions. - /*! @name I/O - */ - //@{ - ///Stream insertion operator - std::ostream& operator<<(std::ostream& os, const vector &v); - ///Stream extraction operator - vector operator>>(std::istream& is, vector& v); - ///Stream insertion operator - std::ostream& operator<<(std::ostream& os, const matrix &M); - ///Stream extraction operator - matrix operator>>(std::istream& is, matrix& v); - //@} +/*! @name I/O + */ +//@{ +///Stream insertion operator +std::ostream& operator<<(std::ostream& os, const vector &v); +///Stream extraction operator +vector operator>>(std::istream& is, vector& v); +///Stream insertion operator +std::ostream& operator<<(std::ostream& os, const matrix &M); +///Stream extraction operator +matrix operator>>(std::istream& is, matrix& v); +//@} - /*! @name Arithmetic - * - * Some arithmetic functions for comfortable syntax. - */ - //@{ +/*! @name Arithmetic + * + * Some arithmetic functions for comfortable syntax. + */ +//@{ - /// Scale a vector. - vector operator*(double a, const vector& v); - /// Euclidean norm of a vector. - double norm(const vector& v); +/// Scale a vector. +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. - matrix inv(const matrix& A); - /// Return copy of transposed matrix. - matrix T(const matrix& A); - /// Trace. - double tr(const matrix& A); - /// Determinant. - double det(matrix& A); - /// L2 condition number, computed with SVD. - double cond(matrix& A); - //@} +/// Scale a matrix. +matrix operator*(double a, const matrix& M); +/// Matrix inverse, computed with LU factorisation. +matrix inv(const matrix& A); +/// Return copy of transposed matrix. +matrix T(const matrix& A); +/// Trace. +double tr(const matrix& A); +/// Determinant. +double det(matrix& A); +/// L2 condition number, computed with SVD. +double cond(matrix& A); +//@} } //Inlined functions namespace kwantix{ - inline double& matrix::operator()(const size_t i, const size_t j){ - try{ - if(LUfactored) - delete LUptr; - if(SVDfactored) - gsl_vector_free(SVD); +inline double& matrix::operator()(const size_t i, const size_t j){ + try{ + if(LUfactored) + delete LUptr; + if(SVDfactored) + gsl_vector_free(SVD); - LUfactored = false; - SVDfactored = false; - return(*gsl_matrix_ptr(A,i-1,j-1)); - } - catch(indexOutOfRange& exc){ - exc.i = i; - exc.j = j; - exc.m = A -> size1; - exc.n = A -> size2; - throw exc; - } + LUfactored = false; + SVDfactored = false; + return(*gsl_matrix_ptr(A,i-1,j-1)); } - inline const double& matrix::operator()(const size_t i, - const size_t j) const{ - try{ - return(*gsl_matrix_const_ptr(A,i-1,j-1)); - } - catch(indexOutOfRange& exc){ - exc.i = i; - exc.j = j; - exc.m = A -> size1; - exc.n = A -> size2; - throw exc; - } - } - - inline double& vector::operator()(const size_t i) { - try{ - return *gsl_vector_ptr(x,i-1); - } - catch(indexOutOfRange& exc){ - exc.i = i; - exc.n = x -> size; - throw exc; - } - } - - inline const double& vector::operator()(const size_t i) const { - try{ - return *gsl_vector_const_ptr(x,i-1); - } - catch(indexOutOfRange& exc){ - exc.i = i; - exc.n = x -> size; - throw exc; - } + catch(indexOutOfRange& exc){ + exc.i = i; + exc.j = j; + exc.m = A -> size1; + exc.n = A -> size2; + throw exc; } } +inline const double& matrix::operator()(const size_t i, + const size_t j) const{ + try{ + return(*gsl_matrix_const_ptr(A,i-1,j-1)); + } + catch(indexOutOfRange& exc){ + exc.i = i; + exc.j = j; + exc.m = A -> size1; + exc.n = A -> size2; + throw exc; + } +} + +inline double& vector::operator()(const size_t i) { + try{ + return *gsl_vector_ptr(x,i-1); + } + catch(indexOutOfRange& exc){ + exc.i = i; + exc.n = x -> size; + throw exc; + } +} + +inline const double& vector::operator()(const size_t i) const { + try{ + return *gsl_vector_const_ptr(x,i-1); + } + catch(indexOutOfRange& exc){ + exc.i = i; + exc.n = x -> size; + throw exc; + } +} +} #endif //__LINALG_HPP__
--- a/src/include/rbf.hpp +++ b/src/include/rbf.hpp @@ -14,223 +14,223 @@ namespace kwantix{ - ///An exception struct for RBFs when attempting to use the wrong dimension. - struct badDimension : public badArgument{ - badDimension() {}; - badDimension(string r, string f, int l) : badArgument(r,f,l){}; - size_t actual_dimension; - size_t given_dimension; - }; +///An exception struct for RBFs when attempting to use the wrong dimension. +struct badDimension : public badArgument{ + badDimension() {}; + badDimension(string r, string f, int l) : badArgument(r,f,l){}; + size_t actual_dimension; + size_t given_dimension; +}; } namespace kwantix{ - ///Base abstract class. - class radial_basis_function : public kwantix::realfunc{ - public: - /// Default constructor, centred at zero by default. - radial_basis_function(); - virtual ~radial_basis_function(); +///Base abstract class. +class radial_basis_function : public kwantix::realfunc{ +public: + /// Default constructor, centred at zero by default. + radial_basis_function(); + virtual ~radial_basis_function(); - ///Specify the dimension of all rbf's. - static void set_dimension(size_t dim); + ///Specify the dimension of all rbf's. + static void set_dimension(size_t dim); - ///Specify the point where this rbf is centred. - radial_basis_function(const point& c); - ///Specify the point where this rbf is centred. - void set_centre(const point& c); + ///Specify the point where this rbf is centred. + radial_basis_function(const point& c); + ///Specify the point where this rbf is centred. + void set_centre(const point& c); - ///Evaluate rbf at x. - double at(const point& x) const; - ///Evaluate rbf at x. - double operator()(const point& x) const; + ///Evaluate rbf at x. + double at(const point& x) const; + ///Evaluate rbf at x. + double operator()(const point& x) const; - ///First partial derivative at x in direction k; - double d(const point& x, size_t k) const; + ///First partial derivative at x in direction k; + double d(const point& x, size_t k) const; - ///Second derivative at x in directions k1 and k2. - double d2(const point& x, size_t k1, size_t k2) const ; + ///Second derivative at x in directions k1 and k2. + double d2(const point& x, size_t k1, size_t k2) const ; - ///An rbf hash function - friend size_t hash_value(const radial_basis_function&); + ///An rbf hash function + friend size_t hash_value(const radial_basis_function&); - //For std::set - bool operator<(const radial_basis_function& phi) const; + //For std::set + bool operator<(const radial_basis_function& phi) const; - protected: +protected: - ///The scalar functions defining each RBF. - virtual double operator()(double r) const = 0; - ///The scalar first derivative defining each RBF. - virtual double d(double r) const = 0; - ///The scalar second derivative defining each RBF. - virtual double d2(double r) const = 0; + ///The scalar functions defining each RBF. + virtual double operator()(double r) const = 0; + ///The scalar first derivative defining each RBF. + virtual double d(double r) const = 0; + ///The scalar second derivative defining each RBF. + virtual double d2(double r) const = 0; - private: +private: - ///In case we attempt to give this rbf an argument of the wrong - ///dimension. - void bad_dimension(string file, int line, size_t dim) const; + ///In case we attempt to give this rbf an argument of the wrong + ///dimension. + void bad_dimension(string file, int line, size_t dim) const; - ///Point on which this rbf is centred. - point centre; + ///Point on which this rbf is centred. + point centre; - ///Dimensionality of all rbfs. - static size_t dimension; - }; + ///Dimensionality of all rbfs. + static size_t dimension; +}; } //Two important subclasses namespace kwantix{ - /// Piecewise smooth RBFs - class piecewise_smooth_rbf : public radial_basis_function{ - public: - piecewise_smooth_rbf(); - piecewise_smooth_rbf(const point &c) :radial_basis_function(c) {;}; - virtual ~piecewise_smooth_rbf(); - }; +/// Piecewise smooth RBFs +class piecewise_smooth_rbf : public radial_basis_function{ +public: + piecewise_smooth_rbf(); + piecewise_smooth_rbf(const point &c) :radial_basis_function(c) {;}; + virtual ~piecewise_smooth_rbf(); +}; - /// \f$C^\infty\f$ RBFs - class c_infty_rbf : public radial_basis_function{ - public: - c_infty_rbf(); - c_infty_rbf(const point &c) :radial_basis_function(c) {;}; +/// \f$C^\infty\f$ RBFs +class c_infty_rbf : public radial_basis_function{ +public: + c_infty_rbf(); + c_infty_rbf(const point &c) :radial_basis_function(c) {;}; - virtual ~c_infty_rbf(); + virtual ~c_infty_rbf(); - static void set_epsilon(double e); - protected: - ///Shape parameter - static double eps; - }; + static void set_epsilon(double e); +protected: + ///Shape parameter + static double eps; +}; } //Specific rbf's follow... namespace kwantix{ - /// \f$r^n\f$ with \f$n\f$ odd - class piecewise_polynomial : public piecewise_smooth_rbf{ - public: - static void set_n(size_t new_n); +/// \f$r^n\f$ with \f$n\f$ odd +class piecewise_polynomial : public piecewise_smooth_rbf{ +public: + static void set_n(size_t new_n); - public: - piecewise_polynomial(){}; - piecewise_polynomial(const point& c) : piecewise_smooth_rbf(c){;}; +public: + piecewise_polynomial(){}; + piecewise_polynomial(const point& c) : piecewise_smooth_rbf(c){;}; - using radial_basis_function::operator(); - using radial_basis_function::d; - using radial_basis_function::d2; + using radial_basis_function::operator(); + using radial_basis_function::d; + using radial_basis_function::d2; - private: - ///Shape parameter - static size_t n; - double operator()(double r) const; - double d(double r) const; - double d2(double r) const; +private: + ///Shape parameter + static size_t n; + double operator()(double r) const; + double d(double r) const; + double d2(double r) const; - }; +}; - ///a common synonym. - typedef piecewise_polynomial conical; +///a common synonym. +typedef piecewise_polynomial conical; } namespace kwantix{ - /// \f$r^n \log(r)\f$ with \f$n\f$ even - class thin_plate_spline : public piecewise_smooth_rbf{ - public: - static void set_n(size_t new_n); +/// \f$r^n \log(r)\f$ with \f$n\f$ even +class thin_plate_spline : public piecewise_smooth_rbf{ +public: + static void set_n(size_t new_n); - public: - thin_plate_spline(){}; - thin_plate_spline(const point& c) : piecewise_smooth_rbf(c){;}; +public: + thin_plate_spline(){}; + thin_plate_spline(const point& c) : piecewise_smooth_rbf(c){;}; - using radial_basis_function::operator(); - using radial_basis_function::d; - using radial_basis_function::d2; + using radial_basis_function::operator(); + using radial_basis_function::d; + using radial_basis_function::d2; - private: - ///Shape parameter - static size_t n; - double operator()(double r) const; - double d(double r) const; - double d2(double r) const; - }; +private: + ///Shape parameter + static size_t n; + double operator()(double r) const; + double d(double r) const; + double d2(double r) const; +}; } namespace kwantix{ - /// \f$\sqrt{1+(\epsilon r)^2}\f$ with \f$\epsilon > 0\f$ - class multiquadric : public c_infty_rbf{ - public: - multiquadric(){}; - multiquadric(const point& c) : c_infty_rbf(c){;}; +/// \f$\sqrt{1+(\epsilon r)^2}\f$ with \f$\epsilon > 0\f$ +class multiquadric : public c_infty_rbf{ +public: + multiquadric(){}; + multiquadric(const point& c) : c_infty_rbf(c){;}; - using radial_basis_function::operator(); - using radial_basis_function::d; - using radial_basis_function::d2; + using radial_basis_function::operator(); + using radial_basis_function::d; + using radial_basis_function::d2; - private: - double operator()(double r) const; - double d(double r) const; - double d2(double r) const; - }; +private: + double operator()(double r) const; + double d(double r) const; + double d2(double r) const; +}; } namespace kwantix{ - /// \f$1/\sqrt{1 + (\epsilon r)^2}\f$ with \f$\epsilon > 0\f$ - class inverse_multiquadric : public c_infty_rbf{ - public: - inverse_multiquadric(){}; - inverse_multiquadric(const point& c) : c_infty_rbf(c){;}; +/// \f$1/\sqrt{1 + (\epsilon r)^2}\f$ with \f$\epsilon > 0\f$ +class inverse_multiquadric : public c_infty_rbf{ +public: + inverse_multiquadric(){}; + inverse_multiquadric(const point& c) : c_infty_rbf(c){;}; - using radial_basis_function::operator(); - using radial_basis_function::d; - using radial_basis_function::d2; + using radial_basis_function::operator(); + using radial_basis_function::d; + using radial_basis_function::d2; - private: - double operator()(double r) const; - double d(double r) const; - double d2(double r) const; - }; +private: + double operator()(double r) const; + double d(double r) const; + double d2(double r) const; +}; } namespace kwantix{ - /// \f$1/(1+(\epsilon r)^2)\f$ with \f$\epsilon > 0 \f$ - class inverse_quadratic : public c_infty_rbf{ - public: - inverse_quadratic(){}; - inverse_quadratic(const point& c) : c_infty_rbf(c){;}; +/// \f$1/(1+(\epsilon r)^2)\f$ with \f$\epsilon > 0 \f$ +class inverse_quadratic : public c_infty_rbf{ +public: + inverse_quadratic(){}; + inverse_quadratic(const point& c) : c_infty_rbf(c){;}; - using radial_basis_function::operator(); - using radial_basis_function::d; - using radial_basis_function::d2; + using radial_basis_function::operator(); + using radial_basis_function::d; + using radial_basis_function::d2; - private: - double operator()(double r) const; - double d(double r) const; - double d2(double r) const; - }; +private: + double operator()(double r) const; + double d(double r) const; + double d2(double r) const; +}; } namespace kwantix{ - /// \f$ e^{- (\epsilon r)^2}\f$ with \f$\epsilon > 0\f$. - class gaussian : public c_infty_rbf{ - public: - gaussian(){}; - gaussian(const point& c) : c_infty_rbf(c){;}; +/// \f$ e^{- (\epsilon r)^2}\f$ with \f$\epsilon > 0\f$. +class gaussian : public c_infty_rbf{ +public: + gaussian(){}; + gaussian(const point& c) : c_infty_rbf(c){;}; - using radial_basis_function::operator(); - using radial_basis_function::d; - using radial_basis_function::d2; + using radial_basis_function::operator(); + using radial_basis_function::d; + using radial_basis_function::d2; - private: - double operator()(double r) const; - double d(double r) const; - double d2(double r) const; - }; +private: + double operator()(double r) const; + double d(double r) const; + double d2(double r) const; +}; }
--- a/src/include/utils.hpp +++ b/src/include/utils.hpp @@ -15,34 +15,34 @@ #include "linalg.hpp" namespace kwantix{ - ///Clears whitespace from front and back of string s. - std::string trim(const std::string& s); +///Clears whitespace from front and back of string s. +std::string trim(const std::string& s); - ///Does map m contain thing? - template<typename K, typename V> - bool contains(const std::map<K,V>& m, K thing); - ///Does set s contain thing? - template<typename E> - bool contains(const std::set<E>& s, E thing); +///Does map m contain thing? +template<typename K, typename V> +bool contains(const std::map<K,V>& m, K thing); +///Does set s contain thing? +template<typename E> +bool contains(const std::set<E>& s, E thing); - ///Does set s1 include set s2? - template<typename E> - bool includes(const std::set<E>& s1, const std::set<E>& s2); +///Does set s1 include set s2? +template<typename E> +bool includes(const std::set<E>& s1, const std::set<E>& s2); - ///Reads matrices from filenames. - kwantix::matrix read_matrix(std::string filename); - ///Reads vectors from filenames. - kwantix::vector read_vector(std::string filename); +///Reads matrices from filenames. +kwantix::matrix read_matrix(std::string filename); +///Reads vectors from filenames. +kwantix::vector read_vector(std::string filename); - /*! \brief Reads map<point,double> from a matrix. - * - *Last column is the value at each point which is represented in - *turn by the rest of the row. - */ - std::map<kwantix::point, double> read_pd_map(std::string filename); +/*! \brief Reads map<point,double> from a matrix. + * + *Last column is the value at each point which is represented in + *turn by the rest of the row. + */ +std::map<kwantix::point, double> read_pd_map(std::string filename); - ///Outputs some information about generic exceptions. - void show_exception(kwantix::error exc); +///Outputs some information about generic exceptions. +void show_exception(kwantix::error exc); }
--- a/src/include/vtkplot.hpp +++ b/src/include/vtkplot.hpp @@ -5,43 +5,40 @@ #include <string> #include <vtkSmartPointer.h> -#include <vtkPolyData.h> -#include "linalg.hpp" +#include "interpolator.hpp" namespace kwantix{ - ///A class for creating VTK plots - class vtkplot{ +///A class for creating VTK plots +template<typename RBF> +class vtkplot{ +public: + /*! \brief Constructor performs a Delaunay triangulation on data + * \param data - An \f$n \times 3\f$ matrix where the first two + * columns are points in the x-y plane on which a + * Delaunay triangulation will be done, and the + * third column is the value at this point. + */ + vtkplot(const interpolator<RBF>& data); + ///Defines the view direction when plotting. + void set_view_direction(const vector& dir); + ///Without changing the triangulation, update the values. + void update_values(const vector& values); + ///Enable offscreen plotting. + void plot_offscreen(bool yesno); + ///Set the base filename used for saving, to which numbers are appended. + void set_base_filename(const std::string& basename); + ///Actually do the plot, whether offscreen or onscreen. + void plot() const; - public: - /*! \brief Constructor performs a Delaunay triangulation on data - * \param data - An \f$n \times 3\f$ matrix where the first two - * columns are points in the x-y plane on which a - * Delaunay triangulation will be done, and the - * third column is the value at this point. - */ - vtkplot(const matrix& data); - ///Defines the view direction when plotting. - void set_view_direction(const vector& dir); - ///Without changing the triangulation, update the values. - void update_values(const vector& values); - ///Enable offscreen plotting. - void plot_offscreen(bool yesno); - ///Set the base filename used for saving, to which numbers are appended. - void set_base_filename(const std::string& basename); - ///Actually do the plot, whether offscreen or onscreen. - void plot() const; - - private: - ///Base filename - string basename; - ///Direction from which to plot - vector view_direction; - ///The plot data - vtkSmartPointer<vtkPolyData> polydata; - ///Whether to plot offscreen or not - bool offscreen; - }; +private: + ///Base filename + string basename; + ///Direction from which to plot + vector view_direction; + ///Whether to plot offscreen or not + bool offscreen; +}; }
--- a/src/interp_values.cpp +++ b/src/interp_values.cpp @@ -4,172 +4,172 @@ namespace kwantix{ - // ************************ 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 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__); +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); - } + 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__); +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); - } + 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__); +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); - } + 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 - { - 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); +} - 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_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 - kwantix::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__); - } - // Else, we're initting with interp_values that are the right size - // and there's no need to do anything, since they're bona fide - // boundary values (parent copy ctor does the work). +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 + kwantix::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__); + } + // Else, we're initting with interp_values that are the right size + // and there's no need to do anything, since they're bona fide + // boundary values (parent copy ctor does the work). +} - //*************** intr_vals stuff ***************************** +//*************** intr_vals stuff ***************************** - intr_values::intr_values(const interp_values& v_in) - : interp_values(v_in) - { - if(v.size() == m+n){ - //Upper part of the vector contains the boundary - kwantix::slice s(1,n); - v = v(s); - } - else if(v.size() != n){ - throw badArgument("Cannot initialise intr_values from given " - "interp_values", __FILE__,__LINE__); - } - // Else, we're initting with interp_values that are the right size - // and there's no need to do anything, since they're bona fide - // boundary values (parent copy ctor does the work). +intr_values::intr_values(const interp_values& v_in) + : interp_values(v_in) +{ + if(v.size() == m+n){ + //Upper part of the vector contains the boundary + kwantix::slice s(1,n); + v = v(s); } + else if(v.size() != n){ + throw badArgument("Cannot initialise intr_values from given " + "interp_values", __FILE__,__LINE__); + } + // Else, we're initting with interp_values that are the right size + // and there's no need to do anything, since they're bona fide + // boundary values (parent copy ctor does the work). +} - //*************** normals stuff ************************************* +//*************** normals stuff ************************************* - normals::normals(size_t rbfs_hash_in, - const map<point, kwantix::vector>& normals_in, - size_t n_in) : rbfs_hash(rbfs_hash_in), n(n_in) - { - map<point, kwantix::vector>::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; +normals::normals(size_t rbfs_hash_in, + const map<point, kwantix::vector>& normals_in, + size_t n_in) : rbfs_hash(rbfs_hash_in), n(n_in) +{ + map<point, kwantix::vector>::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); - } +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; +} - 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 ****************** +//*************** 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) +{ + return v+a; +} - interp_values operator-(double a, const interp_values& v) - { - kwantix::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) +{ + kwantix::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) +{ + return v*a; +} - interp_values operator/(double a, const interp_values& v) - { - kwantix::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) +{ + kwantix::vector w(v.v.size(),a); + return interp_values(w/v.v, v.rbfs_hash,v.n,v.m); +} }
--- a/src/interpolator.cpp +++ b/src/interpolator.cpp @@ -16,739 +16,739 @@ #include <boost/shared_ptr.hpp> namespace kwantix{ - using boost::shared_ptr; +using boost::shared_ptr; + +template<typename RBF> +interpolator<RBF>::interpolator() +{ + initted = false; +} - template<typename RBF> - interpolator<RBF>::interpolator() +template<typename RBF> +interpolator<RBF>::interpolator(shared_ptr<linear_BVP2> bvp) +{ + //Workaround because gdb can't break inside constructors. :-( + init(bvp); +} + +template<typename RBF> +interpolator<RBF>::interpolator(shared_ptr<domain> Omega, + const map<point, double>& Xi) +{ + //Check that Xi matches the domain in size + size_t Omega_size = Omega -> get_interior().size() + + Omega -> get_boundary().size(); + if(Xi.size() != Omega_size) { - initted = false; - } - - template<typename RBF> - interpolator<RBF>::interpolator(shared_ptr<linear_BVP2> bvp) - { - //Workaround because gdb can't break inside constructors. :-( - init(bvp); + badArgument exc; + if(Xi.size() < Omega_size) + exc.reason = "Did not provide enough interpolation data for every " + "point in the given domain."; + else + exc.reason = "Provided more interpolation data than points in the " + "given domain."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } - template<typename RBF> - interpolator<RBF>::interpolator(shared_ptr<domain> Omega, - const map<point, double>& Xi) + //Create a trivial bvp taking into account the given domain. + shared_ptr<Id_op> Id(new Id_op); + shared_ptr<dirichlet_op> D(new dirichlet_op); + map<point, double> f, g; + + //Extract f and g from given information + for(map<point, double>::const_iterator I = Xi.begin(); + I != Xi.end(); I++) { - //Check that Xi matches the domain in size - size_t Omega_size = Omega -> get_interior().size() - + Omega -> get_boundary().size(); - if(Xi.size() != Omega_size) + if(kwantix::contains(Omega -> get_interior(), I -> first)) + f[I -> first] = I -> second; + else if(kwantix::contains(Omega -> get_boundary(), I -> first)) + g[I -> first] = I -> second; + else { badArgument exc; - if(Xi.size() < Omega_size) - exc.reason = "Did not provide enough interpolation data for every " - "point in the given domain."; - else - exc.reason = "Provided more interpolation data than points in the " - "given domain."; + exc.reason = "The interpolation data contains points not in " + "the given domain."; exc.line = __LINE__; exc.file = __FILE__; throw exc; } - - //Create a trivial bvp taking into account the given domain. - shared_ptr<Id_op> Id(new Id_op); - shared_ptr<dirichlet_op> D(new dirichlet_op); - map<point, double> f, g; + } + + shared_ptr<linear_BVP2> new_bvp(new linear_BVP2(Omega, Id, D, f,g)); + + interpolate(new_bvp); +} - //Extract f and g from given information - for(map<point, double>::const_iterator I = Xi.begin(); - I != Xi.end(); I++) - { - if(kwantix::contains(Omega -> get_interior(), I -> first)) - f[I -> first] = I -> second; - else if(kwantix::contains(Omega -> get_boundary(), I -> first)) - g[I -> first] = I -> second; - else - { - badArgument exc; - exc.reason = "The interpolation data contains points not in " - "the given domain."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - } +template<typename RBF> +interpolator<RBF>::interpolator(const map<point, double>& Xi) +{ + interpolate(Xi); +} + +template<typename RBF> +void interpolator<RBF>::interpolate(const map<point, double>& Xi) +{ - shared_ptr<linear_BVP2> new_bvp(new linear_BVP2(Omega, Id, D, f,g)); - - interpolate(new_bvp); - } - - template<typename RBF> - interpolator<RBF>::interpolator(const map<point, double>& Xi) - { - interpolate(Xi); + if(Xi.empty()){//Dude, wtf? + badArgument exc; + exc.reason = "Cannot interpolate if no data is given."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } - template<typename RBF> - void interpolator<RBF>::interpolate(const map<point, double>& Xi) - { - - if(Xi.empty()){//Dude, wtf? - badArgument exc; - exc.reason = "Cannot interpolate if no data is given."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - - //Create a trivial bvp. - shared_ptr<Id_op> Id(new Id_op); - shared_ptr<dirichlet_op> D(new dirichlet_op); - set<point> intr; - set<point> bdry; //empty - map<point, point> nrml; //empty - map<point, double> g; //empty + //Create a trivial bvp. + shared_ptr<Id_op> Id(new Id_op); + shared_ptr<dirichlet_op> D(new dirichlet_op); + set<point> intr; + set<point> bdry; //empty + map<point, point> nrml; //empty + map<point, double> g; //empty - bool dim_set = false; - size_t dimension = 0; - for(map<point,double>::const_iterator i = Xi.begin(); i != Xi.end(); i++){ - if(!dim_set){ - dimension = (i -> first).size(); - dim_set = true; - } - else if(dimension != (i -> first).size()){ - badArgument exc; - exc.reason = "Inconformant dimensions in interpolation data."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - intr.insert( i->first); + bool dim_set = false; + size_t dimension = 0; + for(map<point,double>::const_iterator i = Xi.begin(); i != Xi.end(); i++){ + if(!dim_set){ + dimension = (i -> first).size(); + dim_set = true; } - shared_ptr<domain> Omega(new domain(dimension, intr,bdry,nrml)); - shared_ptr<linear_BVP2> bvp(new linear_BVP2(Omega, Id, D, Xi, g)); - - init(bvp); - } - - template<typename RBF> - void interpolator<RBF>::interpolate(shared_ptr<linear_BVP2> bvp) - { - init(bvp); - } - - template<typename RBF> - void interpolator<RBF>::interpolate(const interp_values& values) - { - if(!initted){ + else if(dimension != (i -> first).size()){ 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.reason = "Inconformant dimensions in interpolation data."; exc.line = __LINE__; exc.file = __FILE__; throw exc; } - coeffs = M.inv(values.v); - - //Precompute the relevant values - recompute_values_vec(); + intr.insert( i->first); } + shared_ptr<domain> Omega(new domain(dimension, intr,bdry,nrml)); + shared_ptr<linear_BVP2> bvp(new linear_BVP2(Omega, Id, D, Xi, g)); + + init(bvp); +} - template<typename RBF> - void interpolator<RBF>::init(shared_ptr<linear_BVP2> bvp) - { - thebvp = bvp; +template<typename RBF> +void interpolator<RBF>::interpolate(shared_ptr<linear_BVP2> bvp) +{ + init(bvp); +} - using std::set; +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 + recompute_values_vec(); +} + +template<typename RBF> +void interpolator<RBF>::init(shared_ptr<linear_BVP2> bvp) +{ + thebvp = bvp; + + using std::set; - shared_ptr<const domain> Omega = bvp -> get_domain(); - set<point> interior = Omega -> get_interior(); - set<point> boundary = Omega -> get_boundary(); - map<point, vector> normals = Omega -> get_normals(); - n = interior.size(); - m = boundary.size(); + shared_ptr<const domain> Omega = bvp -> get_domain(); + set<point> interior = Omega -> get_interior(); + set<point> boundary = Omega -> get_boundary(); + map<point, vector> normals = Omega -> get_normals(); + n = interior.size(); + m = boundary.size(); - vector temp(n+m); - coeffs = temp; - rbfs.reserve(n+m); + vector temp(n+m); + coeffs = temp; + rbfs.reserve(n+m); - RBF::set_dimension(Omega -> get_dimension()); - - set<point>::iterator I; - //Define all the rbfs... - for(I = interior.begin(); I != interior.end(); I++){ - RBF r(*I); - rbfs.push_back(r); - } - for(I = boundary.begin(); I != boundary.end(); I++){ - RBF r(*I); - rbfs.push_back(r); - } + RBF::set_dimension(Omega -> get_dimension()); - //Now define the matrix to be inverted... - matrix Mtemp(n+m,n+m); - M = Mtemp; - shared_ptr<const linear_diff_op2> L = thebvp -> get_linear_diff_op2(); - - shared_ptr<const bdry_diff_op> B = thebvp -> get_bdry_diff_op(); - - I = interior.begin(); - for(size_t i = 1; i <= n; i++){ - for(size_t j = 1; j <= n+m; j++) - M(i,j) = L -> at(rbfs[j-1], *I); - I++; - } + set<point>::iterator I; + //Define all the rbfs... + for(I = interior.begin(); I != interior.end(); I++){ + RBF r(*I); + rbfs.push_back(r); + } + for(I = boundary.begin(); I != boundary.end(); I++){ + RBF r(*I); + rbfs.push_back(r); + } - map<point, vector>::iterator J; - J = normals.begin(); - for(size_t i = n+1; i <= n+m; i++){ - for(size_t j = 1; j <= n+m; j++) - M(i,j) = B -> at(rbfs[j-1], J->first, J->second); - J++; - } - - computecoeffs(); - - ev_precomp = false; - d1_precomp = false; - d2_precomp = false; - - rbfs_hash = hash_value(rbfs); - - kwantix::normals nrmls_tmp(rbfs_hash, normals, n); - nrmls = nrmls_tmp; - - initted = true; - } + //Now define the matrix to be inverted... + matrix Mtemp(n+m,n+m); + M = Mtemp; + shared_ptr<const linear_diff_op2> L = thebvp -> get_linear_diff_op2(); - template<typename RBF> - size_t interpolator<RBF>::hash_value(const std::vector<RBF>& rbfs_in) - { - size_t seed = 0; - - for(size_t i = 0; i < rbfs_in.size(); i++) - boost::hash_combine(seed,rbfs[i]); - - return seed; - } - - // ******************** Evaluations, derivatives ******************** - - // ---------- Precomputation --------------------------- - - template<typename RBF> - void interpolator<RBF>::precompute_ev() - { - if(!initted){ - throw not_initted(__LINE__, __FILE__); - } - - if(!ev_precomp){ - std::vector<size_t> alpha; //empty vector - - //Precompute the RBFs at all points of the domain - matrix phis(n+m,n+m); - set<point>::iterator I; - size_t i; - for(I = thebvp -> get_domain () -> get_interior().begin(), i = 1; - i <= n; i++, I++) - for(size_t j = 0; j < n+m; j++) - phis(i,j+1) = rbfs[j].at(*I); - - for(I = thebvp -> get_domain() -> get_boundary().begin(), i=n+1; - i <= n+m; i++, I++) - for(size_t j = 0; j < n+m; j++) - phis(i,j+1) = rbfs[j].at(*I); - - precomp_rbfs[alpha] = phis; - computecoeffs(); - - ev_precomp = true; - } - } - - template<typename RBF> - void interpolator<RBF>::precompute_d1() - { - if(!initted){ - throw not_initted(__LINE__, __FILE__); - } + shared_ptr<const bdry_diff_op> B = thebvp -> get_bdry_diff_op(); - //Have the RBFs been precomputed yet? - if(!d1_precomp){ - //No, so precompute the RBFs first derivatives at all points of - //the domain. - for(size_t k = 1; k <= thebvp -> get_domain() -> get_dimension(); - k++){ - std::vector<size_t> alpha(k); alpha[k-1]++; - matrix phis(n+m,n+m); - { - set<point>::iterator I; - size_t i; - for(I = thebvp -> get_domain() -> get_interior().begin(), - i = 1; i <= n; i++, I++) - for(size_t j = 0; j < n+m; j++) - phis(i,j+1) = rbfs[j].d(*I,k); - - for(I = thebvp -> get_domain() -> get_boundary().begin(), - i=n+1; i <= n+m; i++, I++) - for(size_t j = 0; j < n+m; j++) - phis(i,j+1) = rbfs[j].d(*I,k); - } - - precomp_rbfs[alpha] = phis; - computecoeffs(); - } - } - d1_precomp = true; + I = interior.begin(); + for(size_t i = 1; i <= n; i++){ + for(size_t j = 1; j <= n+m; j++) + M(i,j) = L -> at(rbfs[j-1], *I); + I++; } + + map<point, vector>::iterator J; + J = normals.begin(); + for(size_t i = n+1; i <= n+m; i++){ + for(size_t j = 1; j <= n+m; j++) + M(i,j) = B -> at(rbfs[j-1], J->first, J->second); + J++; + } + + computecoeffs(); + + ev_precomp = false; + d1_precomp = false; + d2_precomp = false; - template<typename RBF> - void interpolator<RBF>::precompute_d2() - { - if(!initted) - throw not_initted(__LINE__, __FILE__); - if(!d2_precomp){ - size_t dim = thebvp -> get_domain() -> get_dimension(); - //precompute all second derivatives - for(size_t k1 = 1; k1 <= dim; k1++){ - for(size_t k2 = k1; k2 <= dim; k2++){ - std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++; - matrix phis(n+m,n+m); - set<point>::iterator I; - size_t i; - for(I = thebvp -> get_domain () -> get_interior().begin(), i = 1; - i <= n; i++, I++) - for(size_t j = 0; j < n+m; j++) - phis(i,j+1) = rbfs[j].d2(*I,k1,k2); - - for(I = thebvp -> get_domain() -> get_boundary().begin(), i=n+1; - i <= n+m; i++, I++) - for(size_t j = 0; j < n+m; j++) - phis(i,j+1) = rbfs[j].d2(*I,k1,k2); - - precomp_rbfs[alpha] = phis; - computecoeffs(); - } - } - } - d2_precomp = true; - } - - // ----------------- Point evaluation --------------------------- + rbfs_hash = hash_value(rbfs); - template<typename RBF> - double interpolator<RBF>::operator()(const point& p) const - { - return at(p); - } + kwantix::normals nrmls_tmp(rbfs_hash, normals, n); + nrmls = nrmls_tmp; - template<typename RBF> - double interpolator<RBF>::at(const point& p) const - { - if(!initted){ - throw not_initted(__LINE__, __FILE__); - } - std::vector<size_t> alpha; //empty vector + initted = true; +} - //Are we evaluating at a precomputed point in the domain? - if(ev_precomp and - precomp_values[alpha].find(p) != precomp_values[alpha].end()) - return precomp_values[alpha][p]; - - //Else, must compute the value - double result = 0; - for(size_t i = 1; i <= coeffs.size(); i++) - result += coeffs(i)*rbfs[i-1].at(p); - - return result; - } +template<typename RBF> +size_t interpolator<RBF>::hash_value(const std::vector<RBF>& rbfs_in) +{ + size_t seed = 0; - template<typename RBF> - double interpolator<RBF>::d(const point& p, size_t k) const{ - if(!initted){ - throw not_initted(__LINE__, __FILE__); - } - std::vector<size_t> alpha(k); alpha[k-1]++; - - //Are we evaluating at a precomputed point in the domain? - if(d1_precomp and - precomp_values[alpha].find(p) != precomp_values[alpha].end()) - return precomp_values[alpha][p]; - - //Else, must compute the value - double result = 0; - for(size_t i = 1; i <= coeffs.size(); i++) - result += coeffs(i)*rbfs[i-1].d(p,k); - - return result; - } + for(size_t i = 0; i < rbfs_in.size(); i++) + boost::hash_combine(seed,rbfs[i]); - template<typename RBF> - double interpolator<RBF>::d2(const point& p, size_t k1, size_t k2) const{ - if(!initted){ - throw not_initted(__LINE__, __FILE__); - } - std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++; - - //Are we evaluating at a precomputed point in the domain? - if(d2_precomp and - precomp_values[alpha].find(p) != precomp_values[alpha].end()) - return precomp_values[alpha][p]; - - //Else, must compute the value - double result = 0; - for(size_t i = 1; i <= coeffs.size(); i++) - result += coeffs(i)*rbfs[i-1].d2(p,k1,k2); - - return result; + return seed; +} + +// ******************** Evaluations, derivatives ******************** + +// ---------- Precomputation --------------------------- + +template<typename RBF> +void interpolator<RBF>::precompute_ev() +{ + if(!initted){ + throw not_initted(__LINE__, __FILE__); } - // ------------------- Whole domain evaluations ---------------- - - template<typename RBF> - kwantix::badArgument - 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; - return exc; - } - - template<typename RBF> - interp_values interpolator<RBF>::operator()() const - { - return at(); + if(!ev_precomp){ + std::vector<size_t> alpha; //empty vector + + //Precompute the RBFs at all points of the domain + matrix phis(n+m,n+m); + set<point>::iterator I; + size_t i; + for(I = thebvp -> get_domain () -> get_interior().begin(), i = 1; + i <= n; i++, I++) + for(size_t j = 0; j < n+m; j++) + phis(i,j+1) = rbfs[j].at(*I); + + for(I = thebvp -> get_domain() -> get_boundary().begin(), i=n+1; + i <= n+m; i++, I++) + for(size_t j = 0; j < n+m; j++) + phis(i,j+1) = rbfs[j].at(*I); + + precomp_rbfs[alpha] = phis; + computecoeffs(); + + ev_precomp = true; } - - template<typename RBF> - interp_values interpolator<RBF>::at() const - { - if(!ev_precomp) - throw not_precomputed(__LINE__, __FILE__); - - std::vector<size_t> alpha; - return interp_values(precomp_values_vec[alpha],rbfs_hash,n,m); - } - - template<typename RBF> - interp_values interpolator<RBF>::d(size_t k) const - { - if(!d1_precomp) - throw not_precomputed(__LINE__, __FILE__); - - std::vector<size_t> alpha(k); alpha[k-1]++; - return interp_values(precomp_values_vec[alpha],rbfs_hash,n,m); - } - - template<typename RBF> - interp_values interpolator<RBF>::d2(size_t k1, size_t k2) const - { - if(!d2_precomp) - 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,m); - } - - //***************** linear operations *************************** +} - template<typename RBF> - interpolator<RBF> interpolator<RBF>:: - operator+(const interpolator<RBF>& u) const - { - if(this -> rbfs_hash != u.rbfs_hash){ - badArgument exc; - exc.reason = - "Cannot add interpolators on different domains (interior and boundary must match)."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - - interpolator<RBF> out = *this; - out.coeffs = (this -> coeffs) + u.coeffs; - out.recompute_values_vec(); - - return out; - } - - template<typename RBF> - interpolator<RBF> interpolator<RBF>:: - operator-(const interpolator<RBF>& u) const - { - if(this -> rbfs_hash != u.rbfs_hash){ - badArgument exc; - exc.reason = - "Cannot subtract interpolators on different domains (interior and boundary must match)."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - - interpolator<RBF> out = *this; - out.coeffs = (this -> coeffs) - u.coeffs; - out.recompute_values_vec(); - - return out; - } - - template<typename RBF> - interpolator<RBF> interpolator<RBF>::operator*(double a) const - { - interpolator<RBF> u = *this; - u.coeffs = (this -> coeffs)*a; - u.recompute_values_vec(); - return u; - } - - template<typename RBF> - interpolator<RBF> interpolator<RBF>::operator/(double a) const - { - interpolator<RBF> u = *this; - u.coeffs = (this -> coeffs)*(1/a); - u.recompute_values_vec(); - return u; +template<typename RBF> +void interpolator<RBF>::precompute_d1() +{ + if(!initted){ + throw not_initted(__LINE__, __FILE__); } - template<typename RBF> - normals interpolator<RBF>::get_normals() const{ - if(!initted){ - throw badArgument("Interpolator not initialised, so there are " - "no normals!",__FILE__,__LINE__); - } - return nrmls; - } - - //******************* data manipulation ************************** + //Have the RBFs been precomputed yet? + if(!d1_precomp){ + //No, so precompute the RBFs first derivatives at all points of + //the domain. + for(size_t k = 1; k <= thebvp -> get_domain() -> get_dimension(); + k++){ + std::vector<size_t> alpha(k); alpha[k-1]++; + matrix phis(n+m,n+m); + { + set<point>::iterator I; + size_t i; + for(I = thebvp -> get_domain() -> get_interior().begin(), + i = 1; i <= n; i++, I++) + for(size_t j = 0; j < n+m; j++) + phis(i,j+1) = rbfs[j].d(*I,k); + + for(I = thebvp -> get_domain() -> get_boundary().begin(), + i=n+1; i <= n+m; i++, I++) + for(size_t j = 0; j < n+m; j++) + phis(i,j+1) = rbfs[j].d(*I,k); + } - template<typename RBF> - void interpolator<RBF>::set_f(const intr_values& f){ - if(!initted){ - throw not_initted(__LINE__, __FILE__); - } - - std::vector<size_t> alpha; - slice s(1,n); - precomp_values_vec[alpha](s) = f.v; - //precomp_values[alpha] = valsvec2map(precomp_values_vec[alpha]); - - computecoeffs(true); - } - - template<typename RBF> - void interpolator<RBF>::set_g(const bdry_values& g){ - if(!initted){ - throw not_initted(__LINE__, __FILE__); + precomp_rbfs[alpha] = phis; + computecoeffs(); } - - std::vector<size_t> alpha; - slice s(n+1,n+m); - precomp_values_vec[alpha](s) = g.v; - //precomp_values[alpha] = valsvec2map(precomp_values_vec[alpha]); - - computecoeffs(true); - } - - template<typename RBF> - void interpolator<RBF>::set_f(const realfunc& f){ - if(!initted){ - throw not_initted(__LINE__, __FILE__); - } - thebvp -> set_f(f); - computecoeffs(); - } - - template<typename RBF> - void interpolator<RBF>::set_g(const realfunc& g) - { - if(!initted){ - throw not_initted(__LINE__, __FILE__); - } - thebvp -> set_g(g); - computecoeffs(); - } - - template<typename RBF> - void interpolator<RBF>::set_f(const map<point, double>& f) - { - if(!initted){ - throw not_initted(__LINE__, __FILE__); - } - thebvp -> set_f(f); - computecoeffs(); - } - - template<typename RBF> - void interpolator<RBF>::set_g(const map<point, double>& g) - { - if(!initted){ - throw not_initted(__LINE__, __FILE__); - } - thebvp -> set_g(g); - computecoeffs(); } - - 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(); - kwantix::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); - - recompute_values_vec(); - } - - template<typename RBF> - void interpolator<RBF>::recompute_values_vec(){ - 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(); - } + d1_precomp = true; +} - template<typename RBF> - kwantix::badArgument - interpolator<RBF>::not_initted(int line, string file) const - { - badArgument exc; - exc.reason = - "Interpolator can't interpolate without initialisation data."; - exc.line = line; - exc.file = file; - return exc; - } - - template<typename RBF> - void interpolator<RBF>::computecoeffs(bool rhs_defined) - { - kwantix::vector rhs(n+m); - - //Compute the RBF coefficients - if(rhs_defined){ - std::vector<size_t> alpha; - rhs = precomp_values_vec[alpha]; - } - else - { - map<point, double>::const_iterator I; - - I = (thebvp -> get_f()).begin(); - for(size_t i = 1; i <= n; i++){ - rhs(i) = I->second; - I++; - } - I = (thebvp -> get_g()).begin(); - for(size_t i = n+1; i <= n+m; i++){ - rhs(i) = I->second; - I++; - } - } - coeffs = M.inv(rhs); - - //Precompute the values for the derivatives and evaluations that - //have already been computed. - { - typename map<std::vector<size_t>, matrix>::const_iterator I; +template<typename RBF> +void interpolator<RBF>::precompute_d2() +{ + if(!initted) + throw not_initted(__LINE__, __FILE__); + if(!d2_precomp){ + size_t dim = thebvp -> get_domain() -> get_dimension(); + //precompute all second derivatives + for(size_t k1 = 1; k1 <= dim; k1++){ + for(size_t k2 = k1; k2 <= dim; k2++){ + std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++; + matrix phis(n+m,n+m); + set<point>::iterator I; + size_t i; + for(I = thebvp -> get_domain () -> get_interior().begin(), i = 1; + i <= n; i++, I++) + for(size_t j = 0; j < n+m; j++) + phis(i,j+1) = rbfs[j].d2(*I,k1,k2); - for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++){ - kwantix::vector vals = (I->second)*coeffs; - precomp_values_vec[I -> first] = vals; - precomp_values[I -> first] = valsvec2map(vals); + for(I = thebvp -> get_domain() -> get_boundary().begin(), i=n+1; + i <= n+m; i++, I++) + for(size_t j = 0; j < n+m; j++) + phis(i,j+1) = rbfs[j].d2(*I,k1,k2); + + precomp_rbfs[alpha] = phis; + computecoeffs(); } } } + d2_precomp = true; +} - template<typename RBF> - map<point, double> interpolator<RBF>:: - valsvec2map(const kwantix::vector& values) const - { - map<point, double> out; - { - set<point>::const_iterator I; - size_t i; - for(I = thebvp -> get_domain() -> get_interior().begin(), i = 1; - i <= n; I++, i++) - out[*I] = values(i); +// ----------------- Point evaluation --------------------------- + +template<typename RBF> +double interpolator<RBF>::operator()(const point& p) const +{ + return at(p); +} + +template<typename RBF> +double interpolator<RBF>::at(const point& p) const +{ + if(!initted){ + throw not_initted(__LINE__, __FILE__); + } + std::vector<size_t> alpha; //empty vector + + //Are we evaluating at a precomputed point in the domain? + if(ev_precomp and + precomp_values[alpha].find(p) != precomp_values[alpha].end()) + return precomp_values[alpha][p]; + + //Else, must compute the value + double result = 0; + for(size_t i = 1; i <= coeffs.size(); i++) + result += coeffs(i)*rbfs[i-1].at(p); + + return result; +} + +template<typename RBF> +double interpolator<RBF>::d(const point& p, size_t k) const{ + if(!initted){ + throw not_initted(__LINE__, __FILE__); + } + std::vector<size_t> alpha(k); alpha[k-1]++; + + //Are we evaluating at a precomputed point in the domain? + if(d1_precomp and + precomp_values[alpha].find(p) != precomp_values[alpha].end()) + return precomp_values[alpha][p]; + + //Else, must compute the value + double result = 0; + for(size_t i = 1; i <= coeffs.size(); i++) + result += coeffs(i)*rbfs[i-1].d(p,k); + + return result; +} + +template<typename RBF> +double interpolator<RBF>::d2(const point& p, size_t k1, size_t k2) const{ + if(!initted){ + throw not_initted(__LINE__, __FILE__); + } + std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++; + + //Are we evaluating at a precomputed point in the domain? + if(d2_precomp and + precomp_values[alpha].find(p) != precomp_values[alpha].end()) + return precomp_values[alpha][p]; + + //Else, must compute the value + double result = 0; + for(size_t i = 1; i <= coeffs.size(); i++) + result += coeffs(i)*rbfs[i-1].d2(p,k1,k2); + + return result; +} + +// ------------------- Whole domain evaluations ---------------- + +template<typename RBF> +kwantix::badArgument +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; + return exc; +} + +template<typename RBF> +interp_values interpolator<RBF>::operator()() const +{ + return at(); +} + +template<typename RBF> +interp_values interpolator<RBF>::at() const +{ + if(!ev_precomp) + throw not_precomputed(__LINE__, __FILE__); + + std::vector<size_t> alpha; + return interp_values(precomp_values_vec[alpha],rbfs_hash,n,m); +} + +template<typename RBF> +interp_values interpolator<RBF>::d(size_t k) const +{ + if(!d1_precomp) + throw not_precomputed(__LINE__, __FILE__); + + std::vector<size_t> alpha(k); alpha[k-1]++; + return interp_values(precomp_values_vec[alpha],rbfs_hash,n,m); +} + +template<typename RBF> +interp_values interpolator<RBF>::d2(size_t k1, size_t k2) const +{ + if(!d2_precomp) + 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,m); +} + +//***************** linear operations *************************** + +template<typename RBF> +interpolator<RBF> interpolator<RBF>:: +operator+(const interpolator<RBF>& u) const +{ + if(this -> rbfs_hash != u.rbfs_hash){ + badArgument exc; + exc.reason = + "Cannot add interpolators on different domains (interior and boundary must match)."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } + + interpolator<RBF> out = *this; + out.coeffs = (this -> coeffs) + u.coeffs; + out.recompute_values_vec(); + + return out; +} + +template<typename RBF> +interpolator<RBF> interpolator<RBF>:: +operator-(const interpolator<RBF>& u) const +{ + if(this -> rbfs_hash != u.rbfs_hash){ + badArgument exc; + exc.reason = + "Cannot subtract interpolators on different domains (interior and boundary must match)."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } + + interpolator<RBF> out = *this; + out.coeffs = (this -> coeffs) - u.coeffs; + out.recompute_values_vec(); + + return out; +} + +template<typename RBF> +interpolator<RBF> interpolator<RBF>::operator*(double a) const +{ + interpolator<RBF> u = *this; + u.coeffs = (this -> coeffs)*a; + u.recompute_values_vec(); + return u; +} + +template<typename RBF> +interpolator<RBF> interpolator<RBF>::operator/(double a) const +{ + interpolator<RBF> u = *this; + u.coeffs = (this -> coeffs)*(1/a); + u.recompute_values_vec(); + 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 nrmls; +} + +//******************* data manipulation ************************** + +template<typename RBF> +void interpolator<RBF>::set_f(const intr_values& f){ + if(!initted){ + throw not_initted(__LINE__, __FILE__); + } + + std::vector<size_t> alpha; + slice s(1,n); + precomp_values_vec[alpha](s) = f.v; + //precomp_values[alpha] = valsvec2map(precomp_values_vec[alpha]); + + computecoeffs(true); +} + +template<typename RBF> +void interpolator<RBF>::set_g(const bdry_values& g){ + if(!initted){ + throw not_initted(__LINE__, __FILE__); + } + + std::vector<size_t> alpha; + slice s(n+1,n+m); + precomp_values_vec[alpha](s) = g.v; + //precomp_values[alpha] = valsvec2map(precomp_values_vec[alpha]); + + computecoeffs(true); +} + +template<typename RBF> +void interpolator<RBF>::set_f(const realfunc& f){ + if(!initted){ + throw not_initted(__LINE__, __FILE__); + } + thebvp -> set_f(f); + computecoeffs(); +} + +template<typename RBF> +void interpolator<RBF>::set_g(const realfunc& g) +{ + if(!initted){ + throw not_initted(__LINE__, __FILE__); + } + thebvp -> set_g(g); + computecoeffs(); +} + +template<typename RBF> +void interpolator<RBF>::set_f(const map<point, double>& f) +{ + if(!initted){ + throw not_initted(__LINE__, __FILE__); + } + thebvp -> set_f(f); + computecoeffs(); +} + +template<typename RBF> +void interpolator<RBF>::set_g(const map<point, double>& g) +{ + if(!initted){ + throw not_initted(__LINE__, __FILE__); + } + thebvp -> set_g(g); + computecoeffs(); +} + +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(); + kwantix::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); + + recompute_values_vec(); +} + +template<typename RBF> +void interpolator<RBF>::recompute_values_vec(){ + 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; - for(I = thebvp -> get_domain() -> get_boundary().begin(), i = n+1; - i <= n+m; I++, i++) - out[*I] = values(i); - } - return out; + //Pointwise map evaluations we will not compute again. + precomp_values.clear(); +} + +template<typename RBF> +kwantix::badArgument +interpolator<RBF>::not_initted(int line, string file) const +{ + badArgument exc; + exc.reason = + "Interpolator can't interpolate without initialisation data."; + exc.line = line; + exc.file = file; + return exc; +} + +template<typename RBF> +void interpolator<RBF>::computecoeffs(bool rhs_defined) +{ + kwantix::vector rhs(n+m); + + //Compute the RBF coefficients + if(rhs_defined){ + std::vector<size_t> alpha; + rhs = precomp_values_vec[alpha]; } - - //****************** I/O ************************************* - - template<typename RBF> - void interpolator<RBF>::into_os(std::ostream& os) const + else { - kwantix::vector values; - size_t i; - set<point>::const_iterator I; + map<point, double>::const_iterator I; + + I = (thebvp -> get_f()).begin(); + for(size_t i = 1; i <= n; i++){ + rhs(i) = I->second; + I++; + } + I = (thebvp -> get_g()).begin(); + for(size_t i = n+1; i <= n+m; i++){ + rhs(i) = I->second; + I++; + } + } + coeffs = M.inv(rhs); - //If we already have precomputed the values... - std::vector<size_t> alpha; - if(precomp_values_vec.count(alpha)){ - values = precomp_values_vec[alpha]; + //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++){ + kwantix::vector vals = (I->second)*coeffs; + precomp_values_vec[I -> first] = vals; + precomp_values[I -> first] = valsvec2map(vals); } - else{ - kwantix::vector tmp(n+m); - for(i = 1, I=thebvp->get_domain()->get_interior().begin() ; - i<=n; i++, I++) - { - tmp(i) = at(*I); - } + } +} - for(i = n+1, I=thebvp->get_domain()->get_boundary().begin() ; - i<=n+m; i++, I++) - { - tmp(i) = at(*I); - } - - values = tmp; - } +template<typename RBF> +map<point, double> interpolator<RBF>:: +valsvec2map(const kwantix::vector& values) const +{ + map<point, double> out; + { + set<point>::const_iterator I; + size_t i; + for(I = thebvp -> get_domain() -> get_interior().begin(), i = 1; + i <= n; I++, i++) + out[*I] = values(i); + + for(I = thebvp -> get_domain() -> get_boundary().begin(), i = n+1; + i <= n+m; I++, i++) + out[*I] = values(i); + } + return out; +} + +//****************** I/O ************************************* + +template<typename RBF> +void interpolator<RBF>::into_os(std::ostream& os) const +{ + kwantix::vector values; + size_t i; + set<point>::const_iterator I; + //If we already have precomputed the values... + std::vector<size_t> alpha; + if(precomp_values_vec.count(alpha)){ + values = precomp_values_vec[alpha]; + } + else{ + kwantix::vector tmp(n+m); for(i = 1, I=thebvp->get_domain()->get_interior().begin() ; i<=n; i++, I++) { - os << *I << " " << values(i) << std::endl; + tmp(i) = at(*I); } for(i = n+1, I=thebvp->get_domain()->get_boundary().begin() ; i<=n+m; i++, I++) { - os << *I << " " << values(i) << std::endl; + tmp(i) = at(*I); } - + + values = tmp; + } + + for(i = 1, I=thebvp->get_domain()->get_interior().begin() ; + i<=n; i++, I++) + { + os << *I << " " << values(i) << std::endl; } - //Instantiations - template class interpolator<piecewise_polynomial>; - template class interpolator<thin_plate_spline>; - template class interpolator<multiquadric>; - template class interpolator<inverse_multiquadric>; - template class interpolator<inverse_quadratic>; - template class interpolator<gaussian>; + for(i = n+1, I=thebvp->get_domain()->get_boundary().begin() ; + i<=n+m; i++, I++) + { + os << *I << " " << values(i) << std::endl; + } + +} + +//Instantiations +template class interpolator<piecewise_polynomial>; +template class interpolator<thin_plate_spline>; +template class interpolator<multiquadric>; +template class interpolator<inverse_multiquadric>; +template class interpolator<inverse_quadratic>; +template class interpolator<gaussian>; }
--- a/src/linalg.cpp +++ b/src/linalg.cpp @@ -13,405 +13,405 @@ //Matrix stuff namespace kwantix{ - const double matrix::eps = std::numeric_limits<double>::epsilon(); - const double vector::eps = std::numeric_limits<double>::epsilon(); +const double matrix::eps = std::numeric_limits<double>::epsilon(); +const double vector::eps = std::numeric_limits<double>::epsilon(); - // *************** Matrix allocation stuff *************************** +// *************** Matrix allocation stuff *************************** - matrix::matrix(){ - A = gsl_matrix_calloc(1,1); //Allocate 1x1 matrix zero matrix. - LUfactored = false; - SVDfactored = false; - } +matrix::matrix(){ + A = gsl_matrix_calloc(1,1); //Allocate 1x1 matrix zero matrix. + LUfactored = false; + SVDfactored = false; +} + +matrix::matrix(size_t m, size_t n, double fillvalue){ + A = gsl_matrix_alloc(m,n); + gsl_matrix_set_all(A,fillvalue); + LUfactored = false; + SVDfactored = false; +} - matrix::matrix(size_t m, size_t n, double fillvalue){ - A = gsl_matrix_alloc(m,n); - gsl_matrix_set_all(A,fillvalue); - LUfactored = false; - SVDfactored = false; - } +matrix::matrix(gsl_matrix *M){ + A = M; + LUfactored = false; + SVDfactored = false; +} + +matrix::matrix(const matrix& M){ + A = gsl_matrix_alloc(M.A->size1, M.A->size2); + gsl_matrix_memcpy(A,M.A); - matrix::matrix(gsl_matrix *M){ - A = M; - LUfactored = false; - SVDfactored = false; - } + LUfactored = M.LUfactored; + if(LUfactored) + LUptr = new LUmatrix(*M.LUptr); + + condition_number = M.condition_number; - matrix::matrix(const matrix& M){ + SVDfactored = M.SVDfactored; + if(SVDfactored){ + SVD = gsl_vector_alloc(M.SVD->size); + gsl_vector_memcpy(SVD,M.SVD); + } +} + +matrix matrix::operator=(const matrix& M){ + if(this != &M){ + gsl_matrix_free(A); A = gsl_matrix_alloc(M.A->size1, M.A->size2); gsl_matrix_memcpy(A,M.A); + if(LUfactored) + delete LUptr; + LUfactored = M.LUfactored; - if(LUfactored) + if(LUfactored){ LUptr = new LUmatrix(*M.LUptr); - - condition_number = M.condition_number; + } + + if(SVDfactored) + gsl_vector_free(SVD); SVDfactored = M.SVDfactored; if(SVDfactored){ + condition_number = M.condition_number; SVD = gsl_vector_alloc(M.SVD->size); gsl_vector_memcpy(SVD,M.SVD); } } - - matrix matrix::operator=(const matrix& M){ - if(this != &M){ - gsl_matrix_free(A); - A = gsl_matrix_alloc(M.A->size1, M.A->size2); - gsl_matrix_memcpy(A,M.A); - - if(LUfactored) - delete LUptr; - - LUfactored = M.LUfactored; - if(LUfactored){ - LUptr = new LUmatrix(*M.LUptr); - } - - if(SVDfactored) - gsl_vector_free(SVD); + return *this; +} - SVDfactored = M.SVDfactored; - if(SVDfactored){ - condition_number = M.condition_number; - SVD = gsl_vector_alloc(M.SVD->size); - gsl_vector_memcpy(SVD,M.SVD); - } - } - return *this; - } - - matrix::~matrix(){ - if(A != 0) //Has subclass already deleted this matrix? - gsl_matrix_free(A); - if(LUfactored) - delete LUptr; - if(SVDfactored) - gsl_vector_free(SVD); - } +matrix::~matrix(){ + if(A != 0) //Has subclass already deleted this matrix? + gsl_matrix_free(A); + if(LUfactored) + delete LUptr; + if(SVDfactored) + gsl_vector_free(SVD); +} - // *********** Matrix accessor stuff ***************************** - size_t matrix::precsn = 4; +// *********** Matrix accessor stuff ***************************** +size_t matrix::precsn = 4; - size_t matrix::precision() const{ - return precsn; - } +size_t matrix::precision() const{ + return precsn; +} - void matrix::set_precision(size_t p){ - precsn = p; - } +void matrix::set_precision(size_t p){ + precsn = p; +} - size_t matrix::rows() const{ - return A -> size1; - } +size_t matrix::rows() const{ + return A -> size1; +} - size_t matrix::cols() const{ - return A -> size2; - } +size_t matrix::cols() const{ + return A -> size2; +} - vector_view matrix::operator()(const size_t i, const slice &b){ - vector_view x_sub(A,i,b); +vector_view matrix::operator()(const size_t i, const slice &b){ + vector_view x_sub(A,i,b); - if(LUfactored) - delete LUptr; - if(SVDfactored) - gsl_vector_free(SVD); + if(LUfactored) + delete LUptr; + if(SVDfactored) + gsl_vector_free(SVD); - LUfactored = false; - SVDfactored = false; - return x_sub; - } - const vector_view matrix::operator()(const size_t i, const slice &b) const{ - vector_view x_sub(A,i,b); - return x_sub; - } + LUfactored = false; + SVDfactored = false; + return x_sub; +} +const vector_view matrix::operator()(const size_t i, const slice &b) const{ + vector_view x_sub(A,i,b); + return x_sub; +} - vector_view matrix::operator()(const slice &a, const size_t j){ - vector_view x_sub(A,a,j); +vector_view matrix::operator()(const slice &a, const size_t j){ + vector_view x_sub(A,a,j); - if(LUfactored) - delete LUptr; - if(SVDfactored) - gsl_vector_free(SVD); + if(LUfactored) + delete LUptr; + if(SVDfactored) + gsl_vector_free(SVD); - LUfactored = false; - SVDfactored = false; - return x_sub; - } - const vector_view matrix::operator()(const slice &a, const size_t j) const{ - vector_view x_sub(A,a,j); - return x_sub; - } + LUfactored = false; + SVDfactored = false; + return x_sub; +} +const vector_view matrix::operator()(const slice &a, const size_t j) const{ + vector_view x_sub(A,a,j); + return x_sub; +} - // ******************* Matrix arithmetic stuff ****************** - matrix matrix::operator*(const double a) const{ - matrix Z = *this; - gsl_matrix_scale(Z.A,a); - return Z; - } +// ******************* Matrix arithmetic stuff ****************** +matrix matrix::operator*(const double a) const{ + matrix Z = *this; + gsl_matrix_scale(Z.A,a); + return Z; +} - matrix matrix::operator+(const matrix& N) const{ - matrix Z = *this; - try{ - gsl_matrix_add(Z.A,N.A); - } - catch(inconformantSizes& exc){ - exc.n_A = A->size1; - exc.m_A = A->size2; - exc.n_B = N.A -> size1; - exc.m_B = N.A -> size2; - throw exc; - } - return Z; +matrix matrix::operator+(const matrix& N) const{ + matrix Z = *this; + try{ + gsl_matrix_add(Z.A,N.A); } + catch(inconformantSizes& exc){ + exc.n_A = A->size1; + exc.m_A = A->size2; + exc.n_B = N.A -> size1; + exc.m_B = N.A -> size2; + throw exc; + } + return Z; +} - matrix matrix::operator*(const matrix& N)const{ - matrix Z(A->size1,N.A->size2); - try{ - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, A, N.A, 0, Z.A); - } - catch(inconformantSizes& exc){ - exc.n_A = A->size1; - exc.m_A = A->size2; - exc.n_B = N.A -> size1; - exc.m_B = N.A -> size2; - throw exc; - } - return Z; +matrix matrix::operator*(const matrix& N)const{ + matrix Z(A->size1,N.A->size2); + try{ + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, A, N.A, 0, Z.A); } + catch(inconformantSizes& exc){ + exc.n_A = A->size1; + exc.m_A = A->size2; + exc.n_B = N.A -> size1; + exc.m_B = N.A -> size2; + throw exc; + } + return Z; +} - matrix matrix::operator-(const matrix& N) const{ - matrix Z = *this; - try{ - gsl_matrix_sub(Z.A,N.A); - } - catch(inconformantSizes& exc){ - exc.n_A = A->size1; - exc.m_A = A->size2; - exc.n_B = N.A -> size1; - exc.m_B = N.A -> size2; - throw exc; - } - return Z; +matrix matrix::operator-(const matrix& N) const{ + matrix Z = *this; + try{ + gsl_matrix_sub(Z.A,N.A); } + catch(inconformantSizes& exc){ + exc.n_A = A->size1; + exc.m_A = A->size2; + exc.n_B = N.A -> size1; + exc.m_B = N.A -> size2; + throw exc; + } + return Z; +} - vector matrix::operator*(const vector& v) const{ - vector w(A -> size1); - try{ - gsl_blas_dgemv(CblasNoTrans, 1, A, v.x, 0, w.x); - } - catch(inconformantSizes& exc){ - exc.n_A = A->size1; - exc.m_A = A->size2; - exc.n_B = v.x -> size; - exc.m_B = 1; - throw exc; - } - return w; +vector matrix::operator*(const vector& v) const{ + vector w(A -> size1); + try{ + gsl_blas_dgemv(CblasNoTrans, 1, A, v.x, 0, w.x); } + catch(inconformantSizes& exc){ + exc.n_A = A->size1; + exc.m_A = A->size2; + exc.n_B = v.x -> size; + exc.m_B = 1; + throw exc; + } + return w; +} - // ***************** Other arithmetic functions ****************** +// ***************** Other arithmetic functions ****************** - matrix::LUmatrix* matrix::LU() const{ - if(A -> size1 != A -> size2){ - matrixNotSquare exc; - exc.reason = "Cannot LU-factorise non-square matrices."; - exc.file = __FILE__; - exc.line = __LINE__; - throw exc; - } - - if(!LUfactored){ - LUptr = new LUmatrix(A); - gsl_linalg_LU_decomp(LUptr->matrix_ptr(), LUptr->perm_ptr(), - LUptr->sgn_ptr()); - LUfactored = true; - } - return LUptr; +matrix::LUmatrix* matrix::LU() const{ + if(A -> size1 != A -> size2){ + matrixNotSquare exc; + exc.reason = "Cannot LU-factorise non-square matrices."; + exc.file = __FILE__; + exc.line = __LINE__; + throw exc; } - matrix matrix::inv() const{ - if(A -> size1 != A -> size2){ - matrixNotSquare exc; - exc.reason = "Cannot invert non-square matrices."; - exc.file = __FILE__; - exc.line = __LINE__; - throw exc; - } - - if(!LUfactored){ - LUptr = LU(); - LUfactored = true; - } - - matrix Z(LUptr->matrix_ptr()->size1, LUptr->matrix_ptr() -> size2); - - gsl_linalg_LU_invert(LUptr->matrix_ptr(), LUptr->perm_ptr(), Z.A); - return Z; + if(!LUfactored){ + LUptr = new LUmatrix(A); + gsl_linalg_LU_decomp(LUptr->matrix_ptr(), LUptr->perm_ptr(), + LUptr->sgn_ptr()); + LUfactored = true; } + return LUptr; +} - matrix matrix::T() const{ - matrix Z(A->size2, A->size1); - for(size_t i = 1; i <= A->size1; i++) - for(size_t j = 1; j <= A->size2; j++) - Z(j,i) = gsl_matrix_get(A,i-1,j-1); - return Z; +matrix matrix::inv() const{ + if(A -> size1 != A -> size2){ + matrixNotSquare exc; + exc.reason = "Cannot invert non-square matrices."; + exc.file = __FILE__; + exc.line = __LINE__; + throw exc; + } + + if(!LUfactored){ + LUptr = LU(); + LUfactored = true; } + + matrix Z(LUptr->matrix_ptr()->size1, LUptr->matrix_ptr() -> size2); + + gsl_linalg_LU_invert(LUptr->matrix_ptr(), LUptr->perm_ptr(), Z.A); + return Z; +} + +matrix matrix::T() const{ + matrix Z(A->size2, A->size1); + for(size_t i = 1; i <= A->size1; i++) + for(size_t j = 1; j <= A->size2; j++) + Z(j,i) = gsl_matrix_get(A,i-1,j-1); + return Z; +} - double matrix::tr() const{ - if(A -> size1 != A -> size2){ - matrixNotSquare exc; - exc.reason = "Cannot find trace of non-square matrix."; - exc.file = __FILE__; - exc.line = __LINE__; - throw exc; - } - - double result=0; - for(size_t i = 1; i <= A->size1; i++) - result += gsl_matrix_get(A,i,i); - return result; +double matrix::tr() const{ + if(A -> size1 != A -> size2){ + matrixNotSquare exc; + exc.reason = "Cannot find trace of non-square matrix."; + exc.file = __FILE__; + exc.line = __LINE__; + throw exc; } - double matrix::det() const{ - if(A -> size1 != A -> size2){ - matrixNotSquare exc; - exc.reason = "Cannot find determinant of non-square matrices."; - exc.file = __FILE__; - exc.line = __LINE__; - throw exc; - } + double result=0; + for(size_t i = 1; i <= A->size1; i++) + result += gsl_matrix_get(A,i,i); + return result; +} + +double matrix::det() const{ + if(A -> size1 != A -> size2){ + matrixNotSquare exc; + exc.reason = "Cannot find determinant of non-square matrices."; + exc.file = __FILE__; + exc.line = __LINE__; + throw exc; + } - if( !LUfactored ){ - LUptr = LU(); - LUfactored = true; - } + if( !LUfactored ){ + LUptr = LU(); + LUfactored = true; + } - double out = gsl_linalg_LU_det(LUptr->matrix_ptr(),LUptr->sgn()); - return out; + double out = gsl_linalg_LU_det(LUptr->matrix_ptr(),LUptr->sgn()); + return out; +} + +vector matrix::inv(const vector& w) const{ + if(A -> size1 != A -> size2){ + matrixNotSquare exc; + exc.reason = "Cannot invert non-square matrices."; + exc.file = __FILE__; + exc.line = __LINE__; + throw exc; } - vector matrix::inv(const vector& w) const{ - if(A -> size1 != A -> size2){ - matrixNotSquare exc; - exc.reason = "Cannot invert non-square matrices."; - exc.file = __FILE__; - exc.line = __LINE__; - throw exc; - } - - vector v(w.size()); - if( !LUfactored){ - LUptr = LU(); - LUfactored = true; - } - - gsl_linalg_LU_solve(LUptr->matrix_ptr(), LUptr->perm_ptr(), w.x, v.x); - return v; - } - - double matrix::cond() const{ - if(SVDfactored) - return condition_number; - - SVDfactor(); - condition_number = gsl_vector_get(SVD,0) - /gsl_vector_get(SVD,std::max(A->size1, A->size2)-1); - return condition_number; + vector v(w.size()); + if( !LUfactored){ + LUptr = LU(); + LUfactored = true; } - void matrix::SVDfactor() const{ - if(SVDfactored) - return; + gsl_linalg_LU_solve(LUptr->matrix_ptr(), LUptr->perm_ptr(), w.x, v.x); + return v; +} + +double matrix::cond() const{ + if(SVDfactored) + return condition_number; + + SVDfactor(); + condition_number = gsl_vector_get(SVD,0) + /gsl_vector_get(SVD,std::max(A->size1, A->size2)-1); + return condition_number; +} + +void matrix::SVDfactor() const{ + if(SVDfactored) + return; - gsl_matrix * M; - if(A -> size1 < A -> size2){ //Transpose so that m>=n - M = gsl_matrix_alloc(A -> size2, A -> size1); - gsl_matrix_transpose_memcpy(M,A); - } - else{ - M = gsl_matrix_alloc(A -> size1, A -> size2); - gsl_matrix_memcpy(M,A); - } + gsl_matrix * M; + if(A -> size1 < A -> size2){ //Transpose so that m>=n + M = gsl_matrix_alloc(A -> size2, A -> size1); + gsl_matrix_transpose_memcpy(M,A); + } + else{ + M = gsl_matrix_alloc(A -> size1, A -> size2); + gsl_matrix_memcpy(M,A); + } - SVD = gsl_vector_alloc(M -> size2); - gsl_matrix* V = gsl_matrix_alloc(M -> size2, M -> size2); - gsl_vector* work = gsl_vector_alloc(M -> size2); - gsl_linalg_SV_decomp(M,V,SVD,work); - SVDfactored = true; - gsl_vector_free(work); - gsl_matrix_free(V); - gsl_matrix_free(M); - } + SVD = gsl_vector_alloc(M -> size2); + gsl_matrix* V = gsl_matrix_alloc(M -> size2, M -> size2); + gsl_vector* work = gsl_vector_alloc(M -> size2); + gsl_linalg_SV_decomp(M,V,SVD,work); + SVDfactored = true; + gsl_vector_free(work); + gsl_matrix_free(V); + gsl_matrix_free(M); +} } //***************** LUmatrix stuff ********************************** namespace kwantix{ - matrix::LUmatrix::LUmatrix(gsl_matrix* M){ - A = gsl_matrix_alloc(M -> size1, M -> size2); - gsl_matrix_memcpy(A,M); - p = gsl_permutation_alloc(A->size1); - signum = 0; - } +matrix::LUmatrix::LUmatrix(gsl_matrix* M){ + A = gsl_matrix_alloc(M -> size1, M -> size2); + gsl_matrix_memcpy(A,M); + p = gsl_permutation_alloc(A->size1); + signum = 0; +} - matrix::LUmatrix::LUmatrix(const LUmatrix& M){ +matrix::LUmatrix::LUmatrix(const LUmatrix& M){ + p = gsl_permutation_alloc(M.p->size); + gsl_permutation_memcpy(p,M.p); + A = gsl_matrix_alloc(M.A -> size1, M.A -> size2); + gsl_matrix_memcpy(A,M.A); + signum = M.signum; +} + +matrix::LUmatrix matrix::LUmatrix::operator=(const LUmatrix& M){ + if(this != &M){ + gsl_permutation_free(p); p = gsl_permutation_alloc(M.p->size); gsl_permutation_memcpy(p,M.p); + + gsl_matrix_free(A); A = gsl_matrix_alloc(M.A -> size1, M.A -> size2); gsl_matrix_memcpy(A,M.A); - signum = M.signum; - } - matrix::LUmatrix matrix::LUmatrix::operator=(const LUmatrix& M){ - if(this != &M){ - gsl_permutation_free(p); - p = gsl_permutation_alloc(M.p->size); - gsl_permutation_memcpy(p,M.p); - - gsl_matrix_free(A); - A = gsl_matrix_alloc(M.A -> size1, M.A -> size2); - gsl_matrix_memcpy(A,M.A); - - signum = M.signum; - } - return *this; - } - - gsl_matrix* matrix::LUmatrix::matrix_ptr(){ - return A; + signum = M.signum; } - gsl_permutation* matrix::LUmatrix::perm_ptr(){ - return p; - } + return *this; +} - int* matrix::LUmatrix::sgn_ptr(){ - return &signum; - } - - int matrix::LUmatrix::sgn(){ - return signum; - } +gsl_matrix* matrix::LUmatrix::matrix_ptr(){ + return A; +} +gsl_permutation* matrix::LUmatrix::perm_ptr(){ + return p; +} - matrix::LUmatrix::LUmatrix(){ - A = 0; - p = 0; - signum = 0; - } +int* matrix::LUmatrix::sgn_ptr(){ + return &signum; +} + +int matrix::LUmatrix::sgn(){ + return signum; +} - matrix::LUmatrix::~LUmatrix(){ - gsl_permutation_free(p); - p = 0; - gsl_matrix_free(A); - A = 0; - } +matrix::LUmatrix::LUmatrix(){ + A = 0; + p = 0; + signum = 0; +} + +matrix::LUmatrix::~LUmatrix(){ + gsl_permutation_free(p); + p = 0; + gsl_matrix_free(A); + A = 0; +} } @@ -419,556 +419,556 @@ namespace kwantix{ - // **************** Vector allocation stuff ********************* +// **************** Vector allocation stuff ********************* - vector::vector(){ - x = gsl_vector_calloc(1); //Allocate zero vector of size one. - } +vector::vector(){ + x = gsl_vector_calloc(1); //Allocate zero vector of size one. +} - vector::vector(const size_t n, const double fillvalue){ - x = gsl_vector_alloc(n); - gsl_vector_set_all(x,fillvalue); - } +vector::vector(const size_t n, const double fillvalue){ + x = gsl_vector_alloc(n); + gsl_vector_set_all(x,fillvalue); +} + +vector::vector(gsl_vector *y){ + x = y; +} - vector::vector(gsl_vector *y){ - x = y; - } +vector::vector(const gsl_vector *y){ + x = gsl_vector_alloc(y -> size); + gsl_vector_memcpy(x,y); +} - vector::vector(const gsl_vector *y){ - x = gsl_vector_alloc(y -> size); - gsl_vector_memcpy(x,y); - } - - vector::vector(const vector &y){ - x = gsl_vector_alloc(y.x->size); +vector::vector(const vector &y){ + x = gsl_vector_alloc(y.x->size); + gsl_vector_memcpy(x,y.x); +} + +vector& vector::operator=(const vector &y){ + if(this != &y){ + gsl_vector_free(x); + x = gsl_vector_alloc(y.x -> size); gsl_vector_memcpy(x,y.x); } - - vector& vector::operator=(const vector &y){ - if(this != &y){ - gsl_vector_free(x); - x = gsl_vector_alloc(y.x -> size); - gsl_vector_memcpy(x,y.x); - } - return *this; - } + return *this; +} + +vector::~vector(){ + if(x != 0) //Has subclass vector_view already deleted this vector? + gsl_vector_free(x); +} + +// **************** Vector accessor stuff ********************** +size_t vector::precsn = 4; - vector::~vector(){ - if(x != 0) //Has subclass vector_view already deleted this vector? - gsl_vector_free(x); - } +size_t vector::precision() const{ + return precsn; +} + +void vector::set_precision(size_t p){ + precsn = p; +} + +size_t vector::size() const{ + return x->size; +} - // **************** Vector accessor stuff ********************** - size_t vector::precsn = 4; +vector_view vector::operator()(const slice &s) { + vector_view x_sub(x,s); + return x_sub; +} + +const vector_view vector::operator()(const slice &s) const { + vector_view x_sub(x,s); + return x_sub; +} + +//*********************** Vector arithmetic stuff *************** - size_t vector::precision() const{ - return precsn; - } +vector vector::operator*(const double a) const{ + vector v = *this; + gsl_vector_scale(v.x, a); + return v; +} + +vector vector::operator/(const double a) const{ + vector v = *this; + gsl_vector_scale(v.x, 1.0/a); + return v; +} - void vector::set_precision(size_t p){ - precsn = p; +vector vector::operator+(const vector& w) const{ + if(x -> size != w.x->size){ + inconformantSizes exc; + exc.reason = "Cannot add vectors of different sizes."; + exc.file = __FILE__; + exc.line = __LINE__; + exc.n_A = x->size; + exc.n_B = w.x->size; + throw exc; } + vector u = *this; + gsl_vector_add(u.x,w.x); + return u; +} - size_t vector::size() const{ - return x->size; - } - - vector_view vector::operator()(const slice &s) { - vector_view x_sub(x,s); - return x_sub; +vector vector::operator-(const vector& w) const{ + if(x -> size != w.x->size){ + inconformantSizes exc; + exc.reason = "Cannot subtract vectors of different sizes."; + exc.file = __FILE__; + exc.line = __LINE__; + exc.n_A = x->size; + exc.n_B = w.x->size; + throw exc; } + vector u = *this; + gsl_vector_sub(u.x,w.x); + return u; +} - const vector_view vector::operator()(const slice &s) const { - vector_view x_sub(x,s); - return x_sub; +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 arithmetic stuff *************** +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 double a) const{ - vector v = *this; - gsl_vector_scale(v.x, a); - return v; - } +vector vector::operator*(const matrix& M) const +{ + return T(M)* (*this); +} - vector vector::operator/(const double a) const{ - vector v = *this; - gsl_vector_scale(v.x, 1.0/a); - return v; - } +double vector::norm() const +{ + return gsl_blas_dnrm2(x); +} - vector vector::operator+(const vector& w) const{ - if(x -> size != w.x->size){ - inconformantSizes exc; - exc.reason = "Cannot add vectors of different sizes."; - exc.file = __FILE__; - exc.line = __LINE__; - exc.n_A = x->size; - exc.n_B = w.x->size; - throw exc; - } - vector u = *this; - gsl_vector_add(u.x,w.x); - return u; +double vector::operator%(const vector& w) const +{ + double a; + try{ + gsl_blas_ddot(x, w.x, &a); + } + catch(inconformantSizes& exc){ + exc.reason = "Can't dot product vectors of different sizes."; + exc.file = __FILE__; + exc.line = __LINE__; + exc.n_A = x->size; + exc.n_B = w.x -> size; + throw exc; + } + return a; +} + +//Comparison +bool vector::operator==(const vector& w) const{ + if(x -> size != w.x -> size){ + badArgument exc; + exc.reason = "Cannot compare vectors of different sizes."; + exc.file = __FILE__; + exc.line = __LINE__; + throw exc; } - vector vector::operator-(const vector& w) const{ - if(x -> size != w.x->size){ - inconformantSizes exc; - exc.reason = "Cannot subtract vectors of different sizes."; - exc.file = __FILE__; - exc.line = __LINE__; - exc.n_A = x->size; - exc.n_B = w.x->size; - throw exc; - } - vector u = *this; - gsl_vector_sub(u.x,w.x); - return u; - } + for(size_t i = 0; i < x -> size; i++) + if(gsl_fcmp(gsl_vector_get(x,i), gsl_vector_get(w.x,i), eps) != 0) + return false; + + return true; +} - 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; - } +bool vector::operator<(const vector& w) const{ + if(x -> size < w.x -> size) //Smaller vectors go first in this order. + return true; - 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); + for(size_t i = 0; i < x -> size; i++){ + double L = gsl_vector_get(x,i); + double R = gsl_vector_get(w.x,i); + if(L < R ) + return true; + if(L > R ) + return false; } - - 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); - } - catch(inconformantSizes& exc){ - exc.reason = "Can't dot product vectors of different sizes."; - exc.file = __FILE__; - exc.line = __LINE__; - exc.n_A = x->size; - exc.n_B = w.x -> size; - throw exc; - } - return a; - } - - //Comparison - bool vector::operator==(const vector& w) const{ - if(x -> size != w.x -> size){ - badArgument exc; - exc.reason = "Cannot compare vectors of different sizes."; - exc.file = __FILE__; - exc.line = __LINE__; - throw exc; - } - - for(size_t i = 0; i < x -> size; i++) - if(gsl_fcmp(gsl_vector_get(x,i), gsl_vector_get(w.x,i), eps) != 0) - return false; - - return true; - } - - bool vector::operator<(const vector& w) const{ - if(x -> size < w.x -> size) //Smaller vectors go first in this order. - return true; - - for(size_t i = 0; i < x -> size; i++){ - double L = gsl_vector_get(x,i); - double R = gsl_vector_get(w.x,i); - if(L < R ) - return true; - if(L > R ) - return false; - } - return false; //Then vectors are equal. - } + return false; //Then vectors are equal. +} } //************************* Vector view stuff ******************************* namespace kwantix{ - vector_view::vector_view(){ - x = 0; +vector_view::vector_view(){ + x = 0; +} +vector_view::vector_view(gsl_vector* y, const slice &s):vector(new gsl_vector){ + if(s.end() > y->size - 1){ + indexOutOfRange exc; + exc.file = __FILE__; + exc.line = __LINE__; + exc.reason = "invalid vector slice upper range."; + exc.i = s.end(); + exc.n = y->size; + throw exc; } - vector_view::vector_view(gsl_vector* y, const slice &s):vector(new gsl_vector){ - if(s.end() > y->size - 1){ - indexOutOfRange exc; - exc.file = __FILE__; - exc.line = __LINE__; - exc.reason = "invalid vector slice upper range."; - exc.i = s.end(); - exc.n = y->size; - throw exc; - } - x -> size = (s.end() - s.begin())/s.stride()+1; - x -> data = gsl_vector_ptr(y,s.begin()); - x -> stride = s.stride(); - x -> block = 0; - x -> owner = 0; - } + x -> size = (s.end() - s.begin())/s.stride()+1; + x -> data = gsl_vector_ptr(y,s.begin()); + x -> stride = s.stride(); + x -> block = 0; + x -> owner = 0; +} - vector_view::vector_view(gsl_matrix* A, const slice& a, const size_t j ):vector(new gsl_vector){ - if(a.end() > A -> size1 - 1 or j-1 > A -> size2-1){ - indexOutOfRange exc; - exc.file = __FILE__; - exc.line = __LINE__; - exc.reason = "invalid matrix view range."; - exc.i = a.end(); - exc.j = j; - exc.n = A -> size1; - exc.m = A -> size2; - throw exc; - } - x -> size = (a.end() - a.begin())/a.stride()+1; - x -> data = gsl_matrix_ptr(A, a.begin(), j-1); - x -> stride = a.stride()*(A -> tda); //Note that a longer step is - //necessary here. - x -> block = 0; - x -> owner = 0; +vector_view::vector_view(gsl_matrix* A, const slice& a, const size_t j ):vector(new gsl_vector){ + if(a.end() > A -> size1 - 1 or j-1 > A -> size2-1){ + indexOutOfRange exc; + exc.file = __FILE__; + exc.line = __LINE__; + exc.reason = "invalid matrix view range."; + exc.i = a.end(); + exc.j = j; + exc.n = A -> size1; + exc.m = A -> size2; + throw exc; } + x -> size = (a.end() - a.begin())/a.stride()+1; + x -> data = gsl_matrix_ptr(A, a.begin(), j-1); + x -> stride = a.stride()*(A -> tda); //Note that a longer step is + //necessary here. + x -> block = 0; + x -> owner = 0; +} - vector_view::vector_view(gsl_matrix* A, const size_t i, const slice& b ):vector(new gsl_vector){ - if(b.end() > A -> size2 - 1 or i-1 > A-> size1-1){ - indexOutOfRange exc; - exc.file = __FILE__; - exc.line = __LINE__; - exc.reason = "invalid matrix view range."; - exc.i = i; - exc.j = b.end(); - exc.n = A -> size1; - exc.m = A -> size2; - throw exc; - } - x -> size = (b.end() - b.begin())/b.stride()+1; - x -> data = gsl_matrix_ptr(A, i-1, b.begin()); - x -> stride = b.stride(); - x -> block = 0; - x -> owner = 0; +vector_view::vector_view(gsl_matrix* A, const size_t i, const slice& b ):vector(new gsl_vector){ + if(b.end() > A -> size2 - 1 or i-1 > A-> size1-1){ + indexOutOfRange exc; + exc.file = __FILE__; + exc.line = __LINE__; + exc.reason = "invalid matrix view range."; + exc.i = i; + exc.j = b.end(); + exc.n = A -> size1; + exc.m = A -> size2; + throw exc; } + x -> size = (b.end() - b.begin())/b.stride()+1; + x -> data = gsl_matrix_ptr(A, i-1, b.begin()); + x -> stride = b.stride(); + x -> block = 0; + x -> owner = 0; +} - vector_view& vector_view::operator=(const vector& w){ - if(x->size != w.size()){ - badArgument exc; - exc.reason = "Cannot assign to vector view: incomformant sizes."; - exc.file = __FILE__; - exc.line = __LINE__; - throw exc; - } - gsl_vector_memcpy(x,w.x); - return *this; +vector_view& vector_view::operator=(const vector& w){ + if(x->size != w.size()){ + badArgument exc; + exc.reason = "Cannot assign to vector view: incomformant sizes."; + exc.file = __FILE__; + exc.line = __LINE__; + throw exc; } + gsl_vector_memcpy(x,w.x); + return *this; +} - vector_view& vector_view::operator=(const vector_view& w){ - if(x->size != w.size()){ - badArgument exc; - exc.reason = "Cannot assign to vector view: incomformant sizes."; - exc.file = __FILE__; - exc.line = __LINE__; - throw exc; - } - gsl_vector_memcpy(x,w.x); - return *this; +vector_view& vector_view::operator=(const vector_view& w){ + if(x->size != w.size()){ + badArgument exc; + exc.reason = "Cannot assign to vector view: incomformant sizes."; + exc.file = __FILE__; + exc.line = __LINE__; + throw exc; + } + gsl_vector_memcpy(x,w.x); + return *this; - } +} - vector_view::~vector_view(){ - delete x; - x = 0; - } +vector_view::~vector_view(){ + delete x; + x = 0; +} } //Slice stuff namespace kwantix{ - slice::slice(size_t a, size_t b, size_t k) { - set(a,b,k); - } +slice::slice(size_t a, size_t b, size_t k) { + set(a,b,k); +} - slice slice::operator()(size_t a, size_t b, size_t k) { - return set(a,b,k); - } +slice slice::operator()(size_t a, size_t b, size_t k) { + return set(a,b,k); +} - slice slice::set(size_t a, size_t b, size_t k) { - if(b < a){ - badArgument exc; - 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){ - badArgument exc; - 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; - fin = b-1; - str = k; - return *this; +slice slice::set(size_t a, size_t b, size_t k) { + if(b < a){ + badArgument exc; + 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){ + badArgument exc; + 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; + fin = b-1; + str = k; + return *this; +} } //Non-member functions namespace kwantix{ - // ************* I/O ************************ +// ************* I/O ************************ - std::ostream& operator<<(std::ostream& os, const vector &v){ - os.setf(std::ios::scientific); - os << std::setprecision(int(v.precision()) ); - for(size_t i = 1; i <= v.size(); i++){ - os << " " << std::setw(int(v.precision()+6) ) << v(i) ; +std::ostream& operator<<(std::ostream& os, const vector &v){ + os.setf(std::ios::scientific); + os << std::setprecision(int(v.precision()) ); + for(size_t i = 1; i <= v.size(); i++){ + os << " " << std::setw(int(v.precision()+6) ) << v(i) ; + } + return os; +} + +vector operator>>(std::istream& is, vector& v){ + using namespace std; + string s; + list<double> data; + bool colvector = true; + bool shouldbedone = false; + while(getline(is, s)){ + s = kwantix::trim(s); + if(s[0] == '#' or s.size() == 0) //Blank line or comment character + continue; + + stringstream ss; + ss << s; + + double x; + size_t i = 0; + while(ss >> x){ + if( (i > 0 and colvector == false) or (shouldbedone == true)){ + badArgument exc; + exc.reason = "Cannot read vector: bad format in input"; + exc.file = __FILE__; + exc.line = __LINE__; + throw exc; + } + data.push_back(x); + i++; } - return os; + if(i > 1){ + colvector = false; //So it must be a row vector instead + shouldbedone = true; + } + } + + if(data.size() == 0){ + endOfFile exc; + exc.reason = "Cannot read empty vector from input."; + exc.file = __FILE__; + exc.line = __LINE__; + throw exc; } - vector operator>>(std::istream& is, vector& v){ - using namespace std; - string s; - list<double> data; - bool colvector = true; - bool shouldbedone = false; - while(getline(is, s)){ - s = kwantix::trim(s); - if(s[0] == '#' or s.size() == 0) //Blank line or comment character - continue; + vector w(data.size()); + typedef list<double>::iterator LI; + size_t k = 1; + for(LI i = data.begin(); i != data.end(); i++){ + w(k) = *i; + k++; + } + v = w; + return v; +} + +std::ostream& operator<<(std::ostream& os, const matrix& A){ + os.setf(std::ios::scientific); + os << std::setprecision(int(A.precision()) ); + for(size_t i = 1; i <= A.rows(); i++){ + for(size_t j = 1; j <= A.cols(); j++) + os << " " << std::setw(int(A.precision()+6) ) << A(i,j) << " "; + os << std::endl; + } + return os; +} - stringstream ss; - ss << s; +matrix operator>>(std::istream& is, matrix& A){ + using namespace std; + string line, token; + bool rowset = false; + list<double> data; + + size_t rowsize = 0; + size_t rows = 0; + size_t cols = 0; + while(getline(is, line)){ + line = kwantix::trim(line); + //Blank row or comment character. + if(line[0] == '#' or line.length() == 0) + continue; + + stringstream ss_line; + cols = 0; + ss_line << line; + while(ss_line >> token){ + if(token[0] == '#'){ + break; //Rest of line is comment. + } - double x; - size_t i = 0; - while(ss >> x){ - if( (i > 0 and colvector == false) or (shouldbedone == true)){ - badArgument exc; - exc.reason = "Cannot read vector: bad format in input"; - exc.file = __FILE__; - exc.line = __LINE__; - throw exc; - } + //The following may fail on a C++ implementation that doesn't + //obey IEEE arithmetic (IEC 559). We could check for those, + //but do we really want to compile Octave on C++ + //implementations that don't follow IEEE arithmetic? + else if(token == "NaN"){ + double x = std::numeric_limits<double>::quiet_NaN(); + data.push_back(x); + cols++; + } + else if(token == "Inf"){ + double x = std::numeric_limits<double>::infinity(); + data.push_back(x); + cols++; + } + else if(token == "-Inf"){ + double x = -std::numeric_limits<double>::infinity(); data.push_back(x); - i++; + cols++; } - if(i > 1){ - colvector = false; //So it must be a row vector instead - shouldbedone = true; + else if(token == ","){ + ss_line >> token; } - } + + //This also ignores commas and any other token. I think. If + //there's garbage in the token, I have to see what happens + //here. Do we also need to check for garbage? + else{ + double x; + stringstream ss_token; + ss_token << token; + ss_token >> x; + data.push_back(x); + cols++; + } - if(data.size() == 0){ - endOfFile exc; - exc.reason = "Cannot read empty vector from input."; + } + + //First row gives the number of columns, and all successive rows + //must have the same number of elements. + if(!rowset){ + rowset = true; + rowsize = cols; + } + + if (cols != rowsize){ + badArgument exc; + exc.reason = "Cannot read matrix: bad format in input"; exc.file = __FILE__; exc.line = __LINE__; throw exc; } - - vector w(data.size()); - typedef list<double>::iterator LI; - size_t k = 1; - for(LI i = data.begin(); i != data.end(); i++){ - w(k) = *i; + rows++; + } + if(rows == 0){ + endOfFile exc; + exc.reason = "Cannot read empty matrix from input."; + exc.file = __FILE__; + exc.line = __LINE__; + throw exc; + } + + matrix M(rows,cols); + typedef list<double>::iterator LI; + + LI k = data.begin(); + for(size_t i = 1; i <= rows; i++){ + for(size_t j = 1; j <= cols; j++){ + M(i,j) = *k; k++; } - v = w; - return v; - } - - std::ostream& operator<<(std::ostream& os, const matrix& A){ - os.setf(std::ios::scientific); - os << std::setprecision(int(A.precision()) ); - for(size_t i = 1; i <= A.rows(); i++){ - for(size_t j = 1; j <= A.cols(); j++) - os << " " << std::setw(int(A.precision()+6) ) << A(i,j) << " "; - os << std::endl; - } - return os; } - matrix operator>>(std::istream& is, matrix& A){ - using namespace std; - string line, token; - bool rowset = false; - list<double> data; - - size_t rowsize = 0; - size_t rows = 0; - size_t cols = 0; - while(getline(is, line)){ - line = kwantix::trim(line); - //Blank row or comment character. - if(line[0] == '#' or line.length() == 0) - continue; - - stringstream ss_line; - cols = 0; - ss_line << line; - while(ss_line >> token){ - if(token[0] == '#'){ - break; //Rest of line is comment. - } + A = M; + return A; +} + +// *********** non-member arithmetic functions ******************* - //The following may fail on a C++ implementation that doesn't - //obey IEEE arithmetic (IEC 559). We could check for those, - //but do we really want to compile Octave on C++ - //implementations that don't follow IEEE arithmetic? - else if(token == "NaN"){ - double x = std::numeric_limits<double>::quiet_NaN(); - data.push_back(x); - cols++; - } - else if(token == "Inf"){ - double x = std::numeric_limits<double>::infinity(); - data.push_back(x); - cols++; - } - else if(token == "-Inf"){ - double x = -std::numeric_limits<double>::infinity(); - data.push_back(x); - cols++; - } - else if(token == ","){ - ss_line >> token; - } - - //This also ignores commas and any other token. I think. If - //there's garbage in the token, I have to see what happens - //here. Do we also need to check for garbage? - else{ - double x; - stringstream ss_token; - ss_token << token; - ss_token >> x; - data.push_back(x); - cols++; - } +vector operator*(double a, const vector& v) +{ + return v*a; +} +double norm(const vector& v) +{ + return v.norm(); +} - } - - //First row gives the number of columns, and all successive rows - //must have the same number of elements. - if(!rowset){ - rowset = true; - rowsize = cols; - } - - if (cols != rowsize){ - badArgument exc; - exc.reason = "Cannot read matrix: bad format in input"; - exc.file = __FILE__; - exc.line = __LINE__; - throw exc; - } - rows++; - } - if(rows == 0){ - endOfFile exc; - exc.reason = "Cannot read empty matrix from input."; - exc.file = __FILE__; - exc.line = __LINE__; - throw exc; - } - - matrix M(rows,cols); - typedef list<double>::iterator LI; - - LI k = data.begin(); - for(size_t i = 1; i <= rows; i++){ - for(size_t j = 1; j <= cols; j++){ - M(i,j) = *k; - k++; - } - } +matrix operator*(double a, const matrix& M) +{ + return M*a; +} +matrix inv(const matrix& A) +{ + return A.inv(); +} +matrix T(const matrix& A) +{ + return A.T(); +} +double tr(const matrix& A) +{ + return A.tr(); +} +double det(matrix& A) +{ + return A.det(); +} - A = M; - return A; - } - - // *********** non-member arithmetic functions ******************* - - vector operator*(double a, const vector& v) - { - return v*a; - } - double norm(const vector& v) - { - return v.norm(); - } - - matrix operator*(double a, const matrix& M) - { - return M*a; - } - matrix inv(const matrix& A) - { - return A.inv(); - } - matrix T(const matrix& A) - { - return A.T(); - } - double tr(const matrix& A) - { - return A.tr(); - } - double det(matrix& A) - { - return A.det(); - } - - double cond(matrix& A) - { - return A.cond(); - } +double cond(matrix& A) +{ + return A.cond(); +} }
--- a/src/main-Laplace.cpp +++ b/src/main-Laplace.cpp @@ -56,7 +56,7 @@ //Define the domain from given data shared_ptr<domain> Omega(new domain("data/interior.matrix", "data/boundary.matrix", - "data/normals.matrix")); + "data/normals.matrix")); //Define the interior and boundary operators shared_ptr<Laplacian> L(new Laplacian);
--- a/src/main-profiling.cpp +++ b/src/main-profiling.cpp @@ -19,7 +19,7 @@ using namespace std; using namespace kwantix - point p(2),q(2); + point p(2),q(2); p(1) = P1; p(2) = P2; q(1) = Q1; @@ -31,16 +31,16 @@ // out = phi(q); out = sqrt(1 - + sqrt(1+ - (p(1)-q(1))*(p(1)-q(1)) - + (p(2)-q(2))*(p(2)-q(2)) - ) - ); -// out = sqrt(1 -// + sqrt(1+ -// (P1-Q1)*(P1-Q1) + (P2-Q2)*(P2-Q2) -// ) -// ); + + sqrt(1+ + (p(1)-q(1))*(p(1)-q(1)) + + (p(2)-q(2))*(p(2)-q(2)) + ) + ); + // out = sqrt(1 + // + sqrt(1+ + // (P1-Q1)*(P1-Q1) + (P2-Q2)*(P2-Q2) + // ) + // ); } cout << out << endl; return 0;
--- a/src/main-sw-euler.cpp +++ b/src/main-sw-euler.cpp @@ -46,8 +46,8 @@ class Fgen : public realfunc{ public: void set_interps(const interpolator<RBF>& u0, - const interpolator<RBF>& v0, - const interpolator<RBF>& h0); + 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); @@ -104,8 +104,8 @@ 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) {}; + 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: @@ -123,12 +123,12 @@ template<typename RBF> void save_interp(const shared_ptr<domain> Omega, - const interpolator<RBF>& u, - size_t n, string which_interp); + 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); + const boost::shared_ptr<domain> Omega); //********************** Main ****************************************** int main(){ @@ -141,8 +141,8 @@ map<point, double> h_init = kwantix::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")); + "data/circ_bdry.matrix", + "data/circ_nrml.matrix")); shared_ptr<Id_op> I(new Id_op); shared_ptr<iter_neumann> @@ -213,8 +213,8 @@ template<typename RBF> void Fgen<RBF>::set_interps(const interpolator<RBF>& u0, - const interpolator<RBF>& v0, - const interpolator<RBF>& h0){ + const interpolator<RBF>& v0, + const interpolator<RBF>& h0){ set_u(u0); set_v(v0); set_h(h0); @@ -244,30 +244,30 @@ 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) - ); + 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) - ); + 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) + + u(p)*h.d(p,1) + + h(p)*u.d(p,1) + - v(p)*h.d(p,2) + - h(p)*v.d(p,2) - ); + v(p)*h.d(p,2) + + h(p)*v.d(p,2) + ); } iter_neumann::iter_neumann(size_t u_or_v_in) @@ -283,7 +283,7 @@ } double iter_neumann::at(const realfunc &f, const point &p, - const vector &n) const + const vector &n) const { return n(u_or_v)*f(p); } @@ -298,8 +298,8 @@ template<typename RBF> void save_interp(const shared_ptr<domain> Omega, - const interpolator<RBF>& u, - size_t n, string which_interp) + const interpolator<RBF>& u, + size_t n, string which_interp) { using namespace std; string filename = "results/" + which_interp; @@ -320,7 +320,7 @@ slice s(1,2); matrix M(Omega -> get_interior().size() + - Omega -> get_boundary().size(), 3); + 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++){ @@ -339,7 +339,7 @@ template<typename RBF> kwantix::vector at_bdry(const interpolator<RBF>& u, - const boost::shared_ptr<domain> Omega) + const boost::shared_ptr<domain> Omega) { vector out(Omega -> get_boundary().size()); std::set<point>::const_iterator I = Omega -> get_boundary().begin(); @@ -354,7 +354,7 @@ template<typename RBF> void bdry_iter(interpolator<RBF> &u, interpolator<RBF> &v, - const boost::shared_ptr<domain> Omega) + const boost::shared_ptr<domain> Omega) { Gu_or_v<RBF_TYPE> gu(2, Omega), gv(1, Omega); double err = 1;
--- a/src/main-sw-rk4.cpp +++ b/src/main-sw-rk4.cpp @@ -41,8 +41,8 @@ class Fgen { public: void set_interps(const interpolator<RBF>& u0, - const interpolator<RBF>& v0, - const interpolator<RBF>& h0); + 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); @@ -93,8 +93,8 @@ template<typename RBF> void save_interp(const shared_ptr<domain> Omega, - const interpolator<RBF>& u, - size_t n, string which_interp); + const interpolator<RBF>& u, + size_t n, string which_interp); template<typename RBF> void set_bdry_conds(interpolator<RBF> &u, interpolator<RBF> &v); @@ -112,8 +112,8 @@ map<point, double> h_init = kwantix::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")); + "data/circ_bdry.matrix", + "data/circ_nrml.matrix")); shared_ptr<Id_op> I(new Id_op); shared_ptr<dirichlet_op> D(new dirichlet_op); @@ -243,8 +243,8 @@ template<typename RBF> void Fgen<RBF>::set_interps(const interpolator<RBF>& u0, - const interpolator<RBF>& v0, - const interpolator<RBF>& h0){ + const interpolator<RBF>& v0, + const interpolator<RBF>& h0){ set_u(u0); set_v(v0); set_h(h0); @@ -285,8 +285,8 @@ interp_values Fv<RBF>::at() const { return dt*(u()*v.d(1) - + v()*v.d(2) - + g*h.d(2)); + + v()*v.d(2) + + g*h.d(2)); } template<typename RBF> @@ -297,8 +297,8 @@ template<typename RBF> void save_interp(const shared_ptr<domain> Omega, - const interpolator<RBF>& u, - size_t n, string which_interp) + const interpolator<RBF>& u, + size_t n, string which_interp) { using namespace std; string filename = "results/" + which_interp; @@ -311,16 +311,16 @@ { slice s(1,2); matrix M(Omega -> get_interior().size() + - Omega -> get_boundary().size(), 3); + 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++){ + 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++){ + i != Omega -> get_boundary().end(); i++){ M(k,s) = *i; M(k,3) = u(*i); k++; @@ -336,10 +336,10 @@ size_t i,j; for(r = 0, j=1; j <= N; r += 1/(N-1),j++){ for(th = 0,i=1; i <= M; th += 2*pi/(M-1),i++){ - kwantix::vector p(2); - p(1) = r*cos(th); - p(2) = r*sin(th); - out(j,i) = u(p); + kwantix::vector p(2); + p(1) = r*cos(th); + p(2) = r*sin(th); + out(j,i) = u(p); } }
--- a/src/rbf.cpp +++ b/src/rbf.cpp @@ -12,134 +12,134 @@ //Radial basis function stuff namespace kwantix{ - size_t radial_basis_function::dimension = 0; - double c_infty_rbf::eps = 1; - size_t piecewise_polynomial::n = 0; - size_t thin_plate_spline::n = 0; +size_t radial_basis_function::dimension = 0; +double c_infty_rbf::eps = 1; +size_t piecewise_polynomial::n = 0; +size_t thin_plate_spline::n = 0; - radial_basis_function::radial_basis_function(){ - if(dimension != 0){ - point zero(dimension); //Centred by default at the origin. - centre = zero; - } +radial_basis_function::radial_basis_function(){ + if(dimension != 0){ + point zero(dimension); //Centred by default at the origin. + centre = zero; } +} - radial_basis_function::radial_basis_function(const point& c){ - if(c.size() != dimension) - bad_dimension(__FILE__, __LINE__, c.size()); - else - centre = c; - } +radial_basis_function::radial_basis_function(const point& c){ + if(c.size() != dimension) + bad_dimension(__FILE__, __LINE__, c.size()); + else + centre = c; +} - radial_basis_function::~radial_basis_function(){ - //Nothing to destroy! - } +radial_basis_function::~radial_basis_function(){ + //Nothing to destroy! +} - void radial_basis_function::set_centre(const point& c){ - if(c.size () != dimension) - bad_dimension(__FILE__, __LINE__, c.size()); - else - centre = c; - } +void radial_basis_function::set_centre(const point& c){ + if(c.size () != dimension) + bad_dimension(__FILE__, __LINE__, c.size()); + else + centre = c; +} + +void radial_basis_function::set_dimension(size_t dim){ + dimension = dim; +} - void radial_basis_function::set_dimension(size_t dim){ - dimension = dim; - } +void radial_basis_function::bad_dimension(string file, + int line, size_t dim) const{ + kwantix::badDimension exc; + if(dimension == 0) + exc.reason = + "Vector of wrong dimensionality passed to " + "radial basis function. \n" + "(Did you forget to set the dimension?)"; + else + exc.reason = + "Vector of wrong dimensionality passed to " + "radial basis function."; - void radial_basis_function::bad_dimension(string file, - int line, size_t dim) const{ - kwantix::badDimension exc; - if(dimension == 0) - exc.reason = - "Vector of wrong dimensionality passed to " - "radial basis function. \n" - "(Did you forget to set the dimension?)"; - else - exc.reason = - "Vector of wrong dimensionality passed to " - "radial basis function."; + exc.line = line; + exc.file = file; + exc.given_dimension = dim; + exc.actual_dimension = dimension; + throw exc; +} - exc.line = line; - exc.file = file; - exc.given_dimension = dim; - exc.actual_dimension = dimension; +double radial_basis_function::at(const point& x) const{ + return operator()(x); +} + +double radial_basis_function::operator()(const point& x) const{ + if(x.size() != dimension) + bad_dimension(__FILE__, __LINE__, x.size()); + return operator()( norm(x - centre) ); +} + +double radial_basis_function::d(const point& x, size_t k) const{ + if(x.size() != dimension) + bad_dimension(__FILE__, __LINE__, x.size()); + else if(k < 1 or k > dimension){ + kwantix::badArgument exc; + exc.reason = "Cannot differentiate wrt given index: out of bounds."; + exc.line = __LINE__; + exc.file = __FILE__; throw exc; } + else if(x == centre) + return d(0); - double radial_basis_function::at(const point& x) const{ - return operator()(x); - } - - double radial_basis_function::operator()(const point& x) const{ - if(x.size() != dimension) - bad_dimension(__FILE__, __LINE__, x.size()); - return operator()( norm(x - centre) ); - } - - double radial_basis_function::d(const point& x, size_t k) const{ - if(x.size() != dimension) - bad_dimension(__FILE__, __LINE__, x.size()); - else if(k < 1 or k > dimension){ - kwantix::badArgument exc; - exc.reason = "Cannot differentiate wrt given index: out of bounds."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - else if(x == centre) - return d(0); + double r = norm(x - centre); + //Call virtualised derived class function. + double result = d(r)*(x(k) - centre(k))/r; + return result; +} - double r = norm(x - centre); - //Call virtualised derived class function. - double result = d(r)*(x(k) - centre(k))/r; - return result; +double radial_basis_function::d2(const point& x, + size_t k1, size_t k2) const{ + if(x.size() != dimension) + bad_dimension(__FILE__, __LINE__, x.size()); + else if(k1 < 1 or k1 > dimension or k2 < 1 or k2 > dimension){ + kwantix::badArgument exc; + exc.reason = "Cannot differentiate wrt given indices: out of bounds."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } - - double radial_basis_function::d2(const point& x, - size_t k1, size_t k2) const{ - if(x.size() != dimension) - bad_dimension(__FILE__, __LINE__, x.size()); - else if(k1 < 1 or k1 > dimension or k2 < 1 or k2 > dimension){ - kwantix::badArgument exc; - exc.reason = "Cannot differentiate wrt given indices: out of bounds."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - else if(x == centre and k1 == k2) - return d2(0); - else if(x == centre and k1 != k2) - return 0; - else if(k1 != k2){ - double r = norm(x-centre); - double r2 = r*r; - double top = (x(k1) - centre(k1)) * (x(k2)-centre(k2)); - return - top * d2(r) / r2 - - top * d(r) / (r2*r); - } + else if(x == centre and k1 == k2) + return d2(0); + else if(x == centre and k1 != k2) + return 0; + else if(k1 != k2){ double r = norm(x-centre); double r2 = r*r; - double top = x(k1) - centre(k1); top = top*top; - double result = top*d2(r)/r2 + d(r)/r - top*d(r)/(r2*r); - return result; + double top = (x(k1) - centre(k1)) * (x(k2)-centre(k2)); + return + top * d2(r) / r2 - + top * d(r) / (r2*r); + } + double r = norm(x-centre); + double r2 = r*r; + double top = x(k1) - centre(k1); top = top*top; + double result = top*d2(r)/r2 + d(r)/r - top*d(r)/(r2*r); + return result; - } +} - size_t hash_value(const radial_basis_function& phi) - { - size_t seed = 0; - using namespace boost; - for(size_t i = 1; i <= phi.centre.size(); i++) - hash_combine(seed,phi.centre(i)); - hash_combine(seed,typeid(phi).name()); - return seed; - } +size_t hash_value(const radial_basis_function& phi) +{ + size_t seed = 0; + using namespace boost; + for(size_t i = 1; i <= phi.centre.size(); i++) + hash_combine(seed,phi.centre(i)); + hash_combine(seed,typeid(phi).name()); + return seed; +} - bool radial_basis_function::operator<(const radial_basis_function& phi) const - { - return (this->centre) < phi.centre; - } +bool radial_basis_function::operator<(const radial_basis_function& phi) const +{ + return (this->centre) < phi.centre; +} } @@ -147,289 +147,289 @@ namespace kwantix{ - piecewise_smooth_rbf::piecewise_smooth_rbf(){ - //Nothing to create! - } +piecewise_smooth_rbf::piecewise_smooth_rbf(){ + //Nothing to create! +} - piecewise_smooth_rbf::~piecewise_smooth_rbf(){ - //Nothing to destroy! - } +piecewise_smooth_rbf::~piecewise_smooth_rbf(){ + //Nothing to destroy! +} - c_infty_rbf::c_infty_rbf(){ - //Nothing to create! - } +c_infty_rbf::c_infty_rbf(){ + //Nothing to create! +} - c_infty_rbf::~c_infty_rbf(){ - //Nothing to destroy! - } +c_infty_rbf::~c_infty_rbf(){ + //Nothing to destroy! +} - void piecewise_polynomial::set_n(size_t n_new){ - if(n_new % 2 != 1){ - badArgument exc; - exc.reason = "Cannot assign an even n to a piecewise polynomial RBF."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - piecewise_polynomial::n = n_new; +void piecewise_polynomial::set_n(size_t n_new){ + if(n_new % 2 != 1){ + badArgument exc; + exc.reason = "Cannot assign an even n to a piecewise polynomial RBF."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } + piecewise_polynomial::n = n_new; +} - void thin_plate_spline::set_n(size_t n_new){ - if(n_new % 2 != 0){ - badArgument exc; - exc.reason = "Cannot assign an odd n to a thin-plate spline RBF."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - thin_plate_spline::n = n_new; - } +void thin_plate_spline::set_n(size_t n_new){ + if(n_new % 2 != 0){ + badArgument exc; + exc.reason = "Cannot assign an odd n to a thin-plate spline RBF."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } + thin_plate_spline::n = n_new; +} - void c_infty_rbf::set_epsilon(double e){ - if(e <= 0){ - badArgument exc; - exc.reason = "C-infinity RBFs must have a positive epsilon."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - eps = e; +void c_infty_rbf::set_epsilon(double e){ + if(e <= 0){ + badArgument exc; + exc.reason = "C-infinity RBFs must have a positive epsilon."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } + eps = e; +} } // ***************** Specific RBFs ******************************* //Piecewise polynomial namespace kwantix{ - double piecewise_polynomial::operator()(double r) const{ - if(n == 0){ - badArgument exc; - exc.reason = - "Parameter n not set for piecewise_polynomial. \n" - "Use piecewise_polynomial::set_n before evaluating."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - r = fabs(r); +double piecewise_polynomial::operator()(double r) const{ + if(n == 0){ + badArgument exc; + exc.reason = + "Parameter n not set for piecewise_polynomial. \n" + "Use piecewise_polynomial::set_n before evaluating."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } + r = fabs(r); - if(n == 1) - return r; - if(n == 3) - return gsl_pow_3(r); - if(n == 5) - return gsl_pow_5(r); - if(n == 7) - return gsl_pow_7(r); - if(n == 9) - return gsl_pow_9(r); - return pow(r,double(n)); - } + if(n == 1) + return r; + if(n == 3) + return gsl_pow_3(r); + if(n == 5) + return gsl_pow_5(r); + if(n == 7) + return gsl_pow_7(r); + if(n == 9) + return gsl_pow_9(r); + return pow(r,double(n)); +} - double piecewise_polynomial::d(double r) const { - if(n == 0){ - badArgument exc; - exc.reason = - "Parameter n not set for piecewise_polynomial. \n" - "Use piecewise_polynomial::set_n before evaluating."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - r = fabs(r); - if(n == 1) - return 1; - if(n == 3) - return 3*gsl_pow_2(r); - if(n == 5) - return 5*gsl_pow_4(r); - if(n == 7) - return 7*gsl_pow_6(r); - if(n == 9) - return 9*gsl_pow_8(r); - return double(n)*pow(r,double(n-1)); +double piecewise_polynomial::d(double r) const { + if(n == 0){ + badArgument exc; + exc.reason = + "Parameter n not set for piecewise_polynomial. \n" + "Use piecewise_polynomial::set_n before evaluating."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } + r = fabs(r); + if(n == 1) + return 1; + if(n == 3) + return 3*gsl_pow_2(r); + if(n == 5) + return 5*gsl_pow_4(r); + if(n == 7) + return 7*gsl_pow_6(r); + if(n == 9) + return 9*gsl_pow_8(r); + return double(n)*pow(r,double(n-1)); +} - double piecewise_polynomial::d2(double r) const { - if(n == 0){ - badArgument exc; - exc.reason = - "Parameter n not set for piecewise_polynomial. \n" - "Use piecewise_polynomial::set_n before evaluating."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - r = fabs(r); - if(n == 1) - return 0; - if(n == 3) - return 6*r; - if(n == 5) - return 20*gsl_pow_3(r); - if(n == 7) - return 42*pow(r,5); - if(n == 9) - return 72*gsl_pow_7(r); - return double(n*(n-1))*pow(r,double(n-2)); +double piecewise_polynomial::d2(double r) const { + if(n == 0){ + badArgument exc; + exc.reason = + "Parameter n not set for piecewise_polynomial. \n" + "Use piecewise_polynomial::set_n before evaluating."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } + r = fabs(r); + if(n == 1) + return 0; + if(n == 3) + return 6*r; + if(n == 5) + return 20*gsl_pow_3(r); + if(n == 7) + return 42*pow(r,5); + if(n == 9) + return 72*gsl_pow_7(r); + return double(n*(n-1))*pow(r,double(n-2)); +} } //Thin-plate spline namespace kwantix{ - double thin_plate_spline::operator()(double r) const { - if(n == 0){ - badArgument exc; - exc.reason = - "Parameter n not set for thin_plate_spline. \n" - "Use thin_plate_spline::set_n before evaluating."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - r = fabs(r); - if(r == 0) - return 0; - if(n == 2) - return gsl_pow_2(r)*log(r); - if(n == 4) - return gsl_pow_4(r)*log(r); - if(n == 6) - return gsl_pow_6(r)*log(r); - if(n == 8) - return gsl_pow_8(r)*log(r); +double thin_plate_spline::operator()(double r) const { + if(n == 0){ + badArgument exc; + exc.reason = + "Parameter n not set for thin_plate_spline. \n" + "Use thin_plate_spline::set_n before evaluating."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } + r = fabs(r); + if(r == 0) + return 0; + if(n == 2) + return gsl_pow_2(r)*log(r); + if(n == 4) + return gsl_pow_4(r)*log(r); + if(n == 6) + return gsl_pow_6(r)*log(r); + if(n == 8) + return gsl_pow_8(r)*log(r); - return pow(r,double(n))*log(r); - } + return pow(r,double(n))*log(r); +} - double thin_plate_spline::d(double r) const { - if(n == 0){ - badArgument exc; - exc.reason = - "Parameter n not set for thin_plate_spline. \n" - "Use thin_plate_spline::set_n before evaluating."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - r = fabs(r); - if(r == 0) - return 0; - if(n == 2) - return r*(2*log(r) + 1); - if(n == 4) - return gsl_pow_3(r)*(4*log(r) + 1); - if(n == 6) - return gsl_pow_5(r)*(6*log(r) + 1); - if(n == 8) - return gsl_pow_7(r)*(8*log(r) + 1); +double thin_plate_spline::d(double r) const { + if(n == 0){ + badArgument exc; + exc.reason = + "Parameter n not set for thin_plate_spline. \n" + "Use thin_plate_spline::set_n before evaluating."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } + r = fabs(r); + if(r == 0) + return 0; + if(n == 2) + return r*(2*log(r) + 1); + if(n == 4) + return gsl_pow_3(r)*(4*log(r) + 1); + if(n == 6) + return gsl_pow_5(r)*(6*log(r) + 1); + if(n == 8) + return gsl_pow_7(r)*(8*log(r) + 1); - return pow(r,double(n-1))*(double(n)*log(r) + 1); - } + return pow(r,double(n-1))*(double(n)*log(r) + 1); +} - double thin_plate_spline::d2(double r) const { - if(n == 0){ - badArgument exc; - exc.reason = - "Parameter n not set for thin_plate_spline. \n" - "Use thin_plate_spline::set_n before evaluating."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - if (n == 2 and r == 0){ - badDomain exc; - exc.reason = - "For n = 2, thin-plate splines do no have a derivative at zero."; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - if(r == 0) - return 0; +double thin_plate_spline::d2(double r) const { + if(n == 0){ + badArgument exc; + exc.reason = + "Parameter n not set for thin_plate_spline. \n" + "Use thin_plate_spline::set_n before evaluating."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } + if (n == 2 and r == 0){ + badDomain exc; + exc.reason = + "For n = 2, thin-plate splines do no have a derivative at zero."; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } + if(r == 0) + return 0; - if(n == 2) - return 2*log(r) + 3; - if(n == 4) - return gsl_pow_2(r)*(12*log(r) + 7); - if(n == 6) - return gsl_pow_4(r)*(30*log(r) + 11); - if(n == 8) - return gsl_pow_6(r)*(56*log(r) + 15); + if(n == 2) + return 2*log(r) + 3; + if(n == 4) + return gsl_pow_2(r)*(12*log(r) + 7); + if(n == 6) + return gsl_pow_4(r)*(30*log(r) + 11); + if(n == 8) + return gsl_pow_6(r)*(56*log(r) + 15); - return pow(r,double(n-2))*(double(n*n - n)*log(r) + double(2*n - 1)); - } + return pow(r,double(n-2))*(double(n*n - n)*log(r) + double(2*n - 1)); +} } //Multiquadric namespace kwantix{ - double multiquadric::operator()(double r) const { - return sqrt(1 + gsl_pow_2(eps*r)); - } +double multiquadric::operator()(double r) const { + return sqrt(1 + gsl_pow_2(eps*r)); +} - double multiquadric::d(double r)const{ - return eps*eps*fabs(r)/sqrt(1 + gsl_pow_2(eps*r)); - } +double multiquadric::d(double r)const{ + return eps*eps*fabs(r)/sqrt(1 + gsl_pow_2(eps*r)); +} - double multiquadric::d2(double r)const{ - double temp = sqrt(1 + gsl_pow_2(eps*r)); - return - eps*eps/temp - gsl_pow_4(eps)*r*r/gsl_pow_3(temp); - } +double multiquadric::d2(double r)const{ + double temp = sqrt(1 + gsl_pow_2(eps*r)); + return + eps*eps/temp - gsl_pow_4(eps)*r*r/gsl_pow_3(temp); +} } //Inverse multiquadric namespace kwantix{ - double inverse_multiquadric::operator()(double r) const { - return 1/sqrt(1 + gsl_pow_2(eps*r)); - } +double inverse_multiquadric::operator()(double r) const { + return 1/sqrt(1 + gsl_pow_2(eps*r)); +} - double inverse_multiquadric::d(double r) const { - return -eps*eps*fabs(r)/gsl_pow_3(sqrt(1 + gsl_pow_2(eps*r))); - } +double inverse_multiquadric::d(double r) const { + return -eps*eps*fabs(r)/gsl_pow_3(sqrt(1 + gsl_pow_2(eps*r))); +} - double inverse_multiquadric::d2(double r) const { - double temp = sqrt(1 + gsl_pow_2(eps*r)); - return - - eps*eps/gsl_pow_3(temp) + 3*gsl_pow_4(eps)*r*r/gsl_pow_5(temp); - } +double inverse_multiquadric::d2(double r) const { + double temp = sqrt(1 + gsl_pow_2(eps*r)); + return + - eps*eps/gsl_pow_3(temp) + 3*gsl_pow_4(eps)*r*r/gsl_pow_5(temp); +} } //Inverse quadratic namespace kwantix{ - double inverse_quadratic::operator()(double r) const { - return 1/(1 + gsl_pow_2(eps*r)); - } +double inverse_quadratic::operator()(double r) const { + return 1/(1 + gsl_pow_2(eps*r)); +} - double inverse_quadratic::d(double r)const{ - return -2*eps*eps*fabs(r)/gsl_pow_2(1 + gsl_pow_2(eps*r)); - } +double inverse_quadratic::d(double r)const{ + return -2*eps*eps*fabs(r)/gsl_pow_2(1 + gsl_pow_2(eps*r)); +} - double inverse_quadratic::d2(double r)const{ - double temp = 1 + gsl_pow_2(eps*r); - double temp2 = temp*temp, temp3 = temp2*temp; - return - 8*gsl_pow_2(eps*eps*r)/temp3 - 2*eps*eps/temp2; +double inverse_quadratic::d2(double r)const{ + double temp = 1 + gsl_pow_2(eps*r); + double temp2 = temp*temp, temp3 = temp2*temp; + return + 8*gsl_pow_2(eps*eps*r)/temp3 - 2*eps*eps/temp2; - } +} } //Gaussian namespace kwantix{ - double gaussian::operator()(double r) const { - return exp(-gsl_pow_2(eps*r)); - } +double gaussian::operator()(double r) const { + return exp(-gsl_pow_2(eps*r)); +} - double gaussian::d(double r) const { - return -2*eps*eps*fabs(r) * exp(-gsl_pow_2(eps*r)); - } +double gaussian::d(double r) const { + return -2*eps*eps*fabs(r) * exp(-gsl_pow_2(eps*r)); +} - double gaussian::d2(double r) const { - double temp = exp(-gsl_pow_2(eps*r)); - return eps*eps*(4*gsl_pow_2(eps*r) - 2) * temp; +double gaussian::d2(double r) const { + double temp = exp(-gsl_pow_2(eps*r)); + return eps*eps*(4*gsl_pow_2(eps*r) - 2) * temp; - } } +}
--- a/src/utils.cpp +++ b/src/utils.cpp @@ -8,102 +8,102 @@ #include "include/error.hpp" namespace kwantix{ - std::string trim(const std::string& s){ - if(s.length() == 0) - return s; - std::size_t beg = s.find_first_not_of(" \a\b\f\n\r\t\v"); - std::size_t end = s.find_last_not_of(" \a\b\f\n\r\t\v"); - if(beg == std::string::npos) // No non-spaces - return ""; - return std::string(s, beg, end - beg + 1); - } +std::string trim(const std::string& s){ + if(s.length() == 0) + return s; + std::size_t beg = s.find_first_not_of(" \a\b\f\n\r\t\v"); + std::size_t end = s.find_last_not_of(" \a\b\f\n\r\t\v"); + if(beg == std::string::npos) // No non-spaces + return ""; + return std::string(s, beg, end - beg + 1); +} - template<typename K, typename V> - bool contains(const std::map<K,V>& m, K thing){ - return m.find(thing) != m.end(); - } +template<typename K, typename V> +bool contains(const std::map<K,V>& m, K thing){ + return m.find(thing) != m.end(); +} - template<typename E> - bool contains(const std::set<E>& s, E thing){ - return s.find(thing) != s.end(); - } +template<typename E> +bool contains(const std::set<E>& s, E thing){ + return s.find(thing) != s.end(); +} - template<typename E> - bool includes(const std::set<E>& s1, const std::set<E>& s2){ - return std::includes(s2.begin(), s2.end(), s1.begin(), s1.end()); - } +template<typename E> +bool includes(const std::set<E>& s1, const std::set<E>& s2){ + return std::includes(s2.begin(), s2.end(), s1.begin(), s1.end()); +} - using kwantix::vector; - using kwantix::matrix; +using kwantix::vector; +using kwantix::matrix; - matrix read_matrix(std::string filename){ - std::ifstream ifs(filename.c_str()); - if(!ifs){ - kwantix::badArgument exc; - exc.reason = "Cannot open file "; - exc.reason += filename; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - matrix v; - ifs >> v; - return v; +matrix read_matrix(std::string filename){ + std::ifstream ifs(filename.c_str()); + if(!ifs){ + kwantix::badArgument exc; + exc.reason = "Cannot open file "; + exc.reason += filename; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } + matrix v; + ifs >> v; + return v; +} - vector read_vector(std::string filename){ - std::ifstream ifs(filename.c_str()); - if(!ifs){ - kwantix::badArgument exc; - exc.reason = "Cannot open file "; - exc.reason += filename; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - vector v; - ifs >> v; - return v; +vector read_vector(std::string filename){ + std::ifstream ifs(filename.c_str()); + if(!ifs){ + kwantix::badArgument exc; + exc.reason = "Cannot open file "; + exc.reason += filename; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; } + vector v; + ifs >> v; + return v; +} - std::map<kwantix::point, double> read_pd_map(std::string filename){ - std::ifstream ifs(filename.c_str()); - if(!ifs){ - kwantix::badArgument exc; - exc.reason = "Cannot open file "; - exc.reason += filename; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - matrix M; - ifs >> M; +std::map<kwantix::point, double> read_pd_map(std::string filename){ + std::ifstream ifs(filename.c_str()); + if(!ifs){ + kwantix::badArgument exc; + exc.reason = "Cannot open file "; + exc.reason += filename; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } + matrix M; + ifs >> M; - if(M.cols() < 2){ - kwantix::badArgument exc; - exc.reason = - "Input matrix to read_pd_map is too narrow. \n" - "Need at least two columns in the input matrix"; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } + if(M.cols() < 2){ + kwantix::badArgument exc; + exc.reason = + "Input matrix to read_pd_map is too narrow. \n" + "Need at least two columns in the input matrix"; + exc.line = __LINE__; + exc.file = __FILE__; + throw exc; + } - std::map <kwantix::point, double> result; - size_t m = M.cols(); - kwantix::slice s(1,m-1); - for(size_t i = 1; i <= M.rows(); i++) - result[M(i,s)] = M(i,m); - return result; - } + std::map <kwantix::point, double> result; + size_t m = M.cols(); + kwantix::slice s(1,m-1); + for(size_t i = 1; i <= M.rows(); i++) + result[M(i,s)] = M(i,m); + return result; +} - void show_exception(kwantix::error exc){ - using namespace std; +void show_exception(kwantix::error exc){ + using namespace std; - cerr << exc.file << ": ""Caught an exception!" << endl; - cerr << exc.file << ":" << exc.line << " error: "<< exc.reason << endl; - cerr << "From file " << exc.file << endl; - } + cerr << exc.file << ": ""Caught an exception!" << endl; + cerr << exc.file << ":" << exc.line << " error: "<< exc.reason << endl; + cerr << "From file " << exc.file << endl; +} } @@ -111,18 +111,18 @@ #include "include/ddm.hpp" //Instantiations namespace kwantix{ - using boost::shared_ptr; - template bool contains(const std::set<kwantix::point>&, kwantix::point E); - template bool contains(const std::map<kwantix::point, kwantix::vector>& m, - kwantix::point thing); - template bool includes(const std::set<kwantix::point>& s1, - const std::set<kwantix::point>& s2); - template bool contains(const std::map<kwantix::point, - shared_ptr<const kwantix::overlapping_domain> >&, - kwantix::point ); +using boost::shared_ptr; +template bool contains(const std::set<kwantix::point>&, kwantix::point E); +template bool contains(const std::map<kwantix::point, kwantix::vector>& m, + kwantix::point thing); +template bool includes(const std::set<kwantix::point>& s1, + const std::set<kwantix::point>& s2); +template bool contains(const std::map<kwantix::point, + shared_ptr<const kwantix::overlapping_domain> >&, + kwantix::point ); - template bool contains(const std::set<shared_ptr - <const kwantix::overlapping_domain> >&, - shared_ptr<const kwantix::overlapping_domain> E); +template bool contains(const std::set<shared_ptr + <const kwantix::overlapping_domain> >&, + shared_ptr<const kwantix::overlapping_domain> E); }
--- a/src/vtkplot.cpp +++ b/src/vtkplot.cpp @@ -1,3 +1,5 @@ +#include "include/vtkplot.hpp" + #include <vtkPoints.h> #include <vtkPointData.h> #include <vtkDelaunay2D.h> @@ -42,10 +44,11 @@ #include <sstream> #include <iomanip> -using namespace std; +// *************** move this stuff into a class ************************ unsigned int GridSize = 50; +using namespace std; double F(double x, double y, double t) { @@ -118,8 +121,6 @@ ->GetPointData() -> SetScalars(newScalars); } -// *************** move this stuff into a class ************************ - typedef double (*Func3d)(double, double, double); void InitOffscreen(bool offscreen); @@ -129,11 +130,11 @@ void PlotPoints(vtkSmartPointer<vtkDelaunay2D> delaunay, bool offscreen); -int main () -{ - TriangulateTerrain(&F); - return 0; -} +//int main () +//{ +// TriangulateTerrain(&F); +// return 0; +//} void InitOffscreen(bool offscreen) {