Mercurial > hg > kwantix

changeset 36:e32da6c38b59

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Reformat source a bit, remove tabs, use C++0x's auto
author Jordi Gutiérrez Hermoso <jordigh@gmail.com>
date Sun, 14 Feb 2010 19:27:28 -0600
parents 22f78a6faa3e
children e87cb53283b4
files src/bvp.cpp src/ddm.cpp src/include/bvp.hpp src/include/ddm.hpp src/vtkplot.cpp
diffstat 5 files changed, 339 insertions(+), 354 deletions(-) [+]
line wrap: on
line diff
--- a/src/bvp.cpp
+++ b/src/bvp.cpp
@@ -9,72 +9,85 @@
 #include "include/func.hpp"
 
 
-namespace kwantix{
+namespace kwantix
+{
   using std::pair;
   using std::make_pair;
     
   //******************** domain stuff ************************************
 
-  domain::domain(size_t dimension){
+  domain::domain(size_t dimension)
+  {
     dim = dimension;
   }
   domain::domain(size_t dimension, set<point> intr, 
-	 set<point> bdry, map<point, vector> ns){
+                 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){
+  domain::domain(string intr, string bdry, string ns)
+  {
 
     bool intr_empty, bdry_empty;
 
     matrix intr_m ;
-    try{
+    try
+    {
       intr_m = read_matrix(intr);
       intr_empty = false;
     }
-    catch(endOfFile){
+    catch(endOfFile)
+    {
       intr_empty = true;
     }
       
     matrix bdry_m;
-    try{
+    try
+    {
       bdry_m = read_matrix(bdry);
       bdry_empty = false;
     }
-    catch(endOfFile){
+    catch(endOfFile)
+    {
       bdry_empty = true;
     }
     
     matrix ns_m;
-    try{
+    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;	  
+    catch(endOfFile& exc)
+    {
+      if(!bdry_empty)
+      {
+        exc.reason = "Boundary is not empty, so normals cannot be either.";
+        throw exc;    
       }
     }
 
     dim = intr_m.cols();
-    if( bdry_m.cols() != dim and !bdry_empty ){
+    if( bdry_m.cols() != dim and !bdry_empty )
+    {
       badArgument exc;
       exc.reason = 
-	"Wrong parameters for domain from filename. \n"
-	"Dimension of boundary (columns) must equal dimension of interior.";
+        "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){
+    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.";
+        "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;
@@ -83,22 +96,24 @@
     
     if(!intr_empty)
       for(size_t i = 1; i <= intr_m.rows(); i++)
-	add_to_interior( intr_m(i,s) );
+        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) );
+        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));
+        add_to_normals(ns_m(i,s1), ns_m(i, s2));
     }
   }
 
-  domain::~domain(){
+  domain::~domain()
+  {
     //Nothing!
   }
 
   //This clears any data already in the domain and sets the dimension.
-  void domain::set_dimension(size_t dimension){
+  void domain::set_dimension(size_t dimension)
+  {
     dim = dimension;
     interior.clear();
     boundary.clear();
@@ -106,25 +121,30 @@
   }
     
   //Add information to the domain
-  void domain::add_to_interior(const set<point> &intr){
-    for(set<point>::const_iterator 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;
+  void domain::add_to_interior(const set<point> &intr)
+  {
+    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);
     }
   }
 
-  void domain::add_to_interior(const point &intr){
-    if(intr.size() != dim){
+  void domain::add_to_interior(const point &intr)
+  {
+    if(intr.size() != dim)
+    {
       badArgument exc;
       exc.reason = 
-	"Cannot assign to domain's interior: inconformant dimensions."; 
+        "Cannot assign to domain's interior: inconformant dimensions."; 
       exc.line = __LINE__;
       exc.file = __FILE__;
       throw exc;
@@ -132,25 +152,30 @@
     interior.insert(intr);
   }
   
-  void domain::add_to_boundary(const set<point> &bdry){
-    for(set<point>::const_iterator 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;
+  void domain::add_to_boundary(const set<point> &bdry)
+  {
+    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);
     }
   }
 
-  void domain::add_to_boundary(const point &bdry){
-    if(bdry.size() != dim){
+  void domain::add_to_boundary(const point &bdry)
+  {
+    if(bdry.size() != dim)
+    {
       badArgument exc;
       exc.reason = 
-	"Cannot assign to domain's boundary: inconformant dimensions."; 
+        "Cannot assign to domain's boundary: inconformant dimensions."; 
       exc.line = __LINE__;
       exc.file = __FILE__;
       throw exc;
@@ -158,36 +183,42 @@
     boundary.insert(bdry);
   }
 
-  void domain::add_to_normals(const map<point, vector> &ns){
-    for(map<point, vector>::const_iterator 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;
+  void domain::add_to_normals(const map<point, vector> &ns)
+  {
+    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;
+      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);
     }
   }
 
-  void domain::add_to_normals(const point &bdry, const vector &n){
-    if (!kwantix::contains(boundary, bdry)){
+  void domain::add_to_normals(const point &bdry, const vector &n)
+  {
+    if (!kwantix::contains(boundary, bdry))
+    {
       badArgument exc;
       exc.reason = "Bad normal given: must match a point on the boundary.";
       exc.line = __LINE__;
       exc.file = __FILE__;
       throw exc;
     }
-    if(bdry.size() != dim or n.size() != dim){
+    if(bdry.size() != dim or n.size() != dim)
+    {
       badArgument exc;
       exc.reason = "Bad normal given: inconformant dimensions.";
       exc.line = __LINE__;
@@ -199,21 +230,26 @@
   }
 
   //Read that information
-  size_t domain::get_dimension() const{
+  size_t domain::get_dimension() const
+  {
     return dim;
   }
-  const set<point>& domain::get_interior() const{
+  const set<point>& domain::get_interior() const
+  {
     return interior;
   }
-  const set<point>& domain::get_boundary() const{
+  const set<point>& domain::get_boundary() const
+  {
     return boundary;
   }
-  const map<point, vector>& domain::get_normals() const{
+  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{
+  bool domain::contains(const point& p) const
+  {
     if(kwantix::contains(interior, p))
       return true;
     if(kwantix::contains(boundary, p))
@@ -223,112 +259,81 @@
   
   //**************** BVP stuff *********************************
 
-  //FIXME: Try to templatise this later
-  BVP::BVP(shared_ptr<const domain> O, 
-	   shared_ptr<const diff_op> L_in, 
-	   shared_ptr<const bdry_diff_op> B_in, 
-	   const realfunc& f_in, 
-	   const realfunc& g_in){
-    Omega = O;
-    L = L_in;
-    B = B_in;
-    set_f(f_in);
-    set_g(g_in);
-  }
-  BVP::BVP(shared_ptr<const domain> O, 
-	   shared_ptr<const diff_op> L_in, 
-	   shared_ptr<const bdry_diff_op> B_in, 
-	   const realfunc& f_in, 
-	   const map<point,double>& g_in){
-    Omega = O;
-    L = L_in;
-    B = B_in;
-    set_f(f_in);
-    set_g(g_in);
-  }
-  BVP::BVP(shared_ptr<const domain> O, 
-	   shared_ptr<const diff_op> L_in, 
-	   shared_ptr<const bdry_diff_op> B_in, 
-	   const map<point,double>& f_in, 
-	   const realfunc& g_in){
-    Omega = O;
-    L = L_in;
-    B = B_in;
-    set_f(f_in);
-    set_g(g_in);
-  }
-  BVP::BVP(shared_ptr<const domain> O, 
-	   shared_ptr<const diff_op> L_in, 
-	   shared_ptr<const bdry_diff_op> B_in, 
-	   const map<point,double>& f_in, 
-	   const map<point,double>& g_in){
-    Omega = O;
-    L = L_in;
-    B = B_in;
-    set_f(f_in);
-    set_g(g_in);
-  }
-
-
-  shared_ptr<const domain> BVP::get_domain() const{
+  shared_ptr<const domain> BVP::get_domain() const
+  {
     return Omega;
   }
-  shared_ptr<const diff_op> BVP::get_diff_op() const{
+  shared_ptr<const diff_op> BVP::get_diff_op() const
+  {
     return L;
   }
-  shared_ptr<const bdry_diff_op> BVP::get_bdry_diff_op() const{
+  shared_ptr<const bdry_diff_op> BVP::get_bdry_diff_op() const
+  {
     return B;
   }  
-  const map<point, double>& BVP::get_f() const{
+  const map<point, double>& BVP::get_f() const
+  {
     return f;
   }
-  const map<point, double>& BVP::get_g() const{
+  const map<point, double>& BVP::get_g() const
+  {
     return g;
   }
 
-  void BVP::set_f(const map<point, double>& f_in){
-    for( map<point,double>::const_iterator I = f_in.begin();
-	 I != f_in.end(); I++){
-      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;
+  void BVP::set_f(const map<point, double>& f_in)
+  {
+    for( auto I = f_in.begin(); I != f_in.end(); I++)
+    {
+      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;
     }
   }
-  void BVP::set_g(const map<point, double>& g_in){
-    for( map<point,double>::const_iterator 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;
+
+  
+  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){
-    for(set<point>::iterator I = Omega->get_interior().begin();
-        I != Omega->get_interior().end(); I++)
+  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);
+    }
   }
 
-  void BVP::set_g(const realfunc &g_in){
-    for(set<point>::iterator I = Omega->get_boundary().begin();
-	I != Omega->get_boundary().end(); I++)
+  void BVP::set_g(const realfunc &g_in)
+  {
+    for(auto I = Omega->get_boundary().begin(); I != Omega->get_boundary().end(); I++)
+    {
       g[*I] = g_in(*I);
+    }
   }
 
-  shared_ptr<const linear_diff_op2> linear_BVP2::get_linear_diff_op2() const{
+  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
@@ -12,7 +12,8 @@
 //debug
 #include <iostream>
 
-namespace kwantix{
+namespace kwantix
+{
   using namespace std;
   using kwantix::vector;
   using boost::shared_ptr;
@@ -21,8 +22,8 @@
   //************* 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) 
+                                     shared_ptr<const bdry_diff_op> Bprime_in,
+                                     const set<point>& ibps) 
   {
     B = B_in ;
     intr_bdry_pts = ibps;
@@ -39,7 +40,7 @@
   }
   
   double ddm_bdry_diff_op::at(const realfunc &f, const point &p, 
-			      const vector &n) const
+                              const vector &n) const
   {
     if( kwantix::contains(intr_bdry_pts, p))     
       return Bprime -> at(f,p,n);
@@ -50,14 +51,14 @@
   // ************ ddm stuff ******************************
 
   ddm::ddm(const set<shared_ptr<const domain> >& ds, 
-		shared_ptr<const BVP> thebvp)
+           shared_ptr<const BVP> thebvp)
   {
     //Gotta check this is actually a domain decomposition...
     set<point> union_interior;
     set<point> union_boundary;
 
-    for(set<shared_ptr<const domain> >::iterator i = ds.begin();
-	i != ds.end(); i++){
+    for(auto i = ds.begin(); i != ds.end(); i++)
+    {
       set<point> intr = (*i) -> get_interior();
       set<point> bdry = (*i) -> get_boundary();
       union_interior.insert(intr.begin(), intr.end());
@@ -69,23 +70,25 @@
     set<point> interior = Omega -> get_interior();
     set<point> boundary = Omega -> get_boundary();
 
-    if( interior != union_interior){
+    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.";
+        "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) ){
+    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 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;
@@ -94,7 +97,7 @@
     domains = ds;
     tolerance = 1e-6;
 
- }
+  }
 
   ddm::~ddm()
   {
@@ -115,13 +118,13 @@
 
   template <typename RBF> additive_schwarz_ddm<RBF>::
   additive_schwarz_ddm(const set<shared_ptr<const domain> >& ds, 
-		       const shared_ptr<const linear_BVP2> thebvp) : 
+                       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());
+    const linear_diff_op2>(this -> bvp -> get_diff_op());
 
     shared_ptr<const bdry_diff_op> 
       B = this -> bvp -> get_bdry_diff_op();
@@ -130,45 +133,45 @@
     map<point, double> global_g = this -> bvp -> get_g();
     
     //Define a bvp for each domain and assign it an interpolator. 
-    for(set<shared_ptr<const domain> >::iterator 
-	  i = this -> domains.begin(); i != this -> domains.end(); i++){
+    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);
+        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(set<point>::iterator spi = this_intr.begin(); 
-	  spi != this_intr.end(); spi++)
-	this_f[*spi] = global_f[*spi];
+      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(set<point>::iterator 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. 
-	}
+      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));
+        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) );
+        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));
+        rbf_ptr( new interpolator<RBF>(this_bvp) );
       phis[this_domain] = rbf_ptr;
     }
     
@@ -185,7 +188,7 @@
     //domains and use the interpolator from that domain to evaluate.
    
     set<shared_ptr<const interpolator<RBF> > >
-    relevant_interpolators = which_interps(p);
+      relevant_interpolators = which_interps(p);
     
     if(relevant_interpolators.size() != 0)//p is in one of the
                                           //domains.
@@ -199,21 +202,21 @@
     //Uh, just begin with closest being some point in the whatever
     //domain. 
     point c = *(((*( this->domains.begin() )) 
-		 -> get_interior()).begin() );
+                 -> get_interior()).begin() );
 
     //Search each domain's interior and boundary. Can't do better than
     //exhaustive search. 
-    for(set<shared_ptr<const domain> >::iterator i = this -> domains.begin();
-	i != this -> domains.end(); i++){
+    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;
+          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;
+          j != (*i) -> get_boundary().end(); j++)
+        if( norm(*j - p) < norm(c - p))
+          c = *j;
     }
     
     return avg_interp( which_interps(c), p);
@@ -226,15 +229,15 @@
   {
     set<shared_ptr<const interpolator<RBF> > > relevant_interpolators;
     for(set<shared_ptr<const domain> >::iterator i = this -> domains.begin();
-	i != this -> domains.end(); i++){
+        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.
+        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;
@@ -243,17 +246,17 @@
   template <typename RBF>
   double additive_schwarz_ddm<RBF>::
   avg_interp(set<shared_ptr<const interpolator<RBF> > > relevant_interpolators,
-	     const point& p) const
+             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;
+    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> 
@@ -274,28 +277,28 @@
 
       //Each domain will have new g's.
       map<shared_ptr<const overlapping_domain>, map<point, double> > 
-	new_bdry_assts;
+        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){
+            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))
-		     );
-	  }
+            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]);
+            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);
@@ -311,19 +314,19 @@
   {
     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.
+          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);
-	  }
+      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() );
@@ -337,33 +340,35 @@
   //************** 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)
+    :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)
+                     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 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;
-	}
+          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 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;
   }
@@ -372,7 +377,7 @@
     :domain(dimension){}
   overlapping_domain::
   overlapping_domain(size_t dimension, set<point> intr, 
-		     set<point> bdry, map<point, vector> ns)
+                     set<point> bdry, map<point, vector> ns)
     :domain(dimension, intr, bdry, ns){}
     
     
@@ -385,7 +390,8 @@
   shared_ptr<const overlapping_domain>
   overlapping_domain::which_domain(const point& p) const
   {
-    if(!kwantix::contains(boundary_assignments, p)){
+    if( !kwantix::contains(boundary_assignments, p) )
+    {
       shared_ptr<const overlapping_domain> zero;
       return zero;
     }
@@ -405,20 +411,22 @@
   void 
   set_overlapper_info(set<shared_ptr<overlapping_domain> > domains)
   {
-    for(set<shared_ptr<overlapping_domain> >::iterator d = domains.begin();
-	d != domains.end(); d++){
-      for(set<point>::iterator p = (*d) -> get_boundary().begin();
-	  p != (*d) -> get_boundary().end(); p++){
-	for(set<shared_ptr<overlapping_domain> >::iterator 
-	      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.
-	  }
+    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++)
+        {
+          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.
+          }
+        }
       }
     }
   }
--- a/src/include/bvp.hpp
+++ b/src/include/bvp.hpp
@@ -36,7 +36,7 @@
     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);
+           set<point> bdry, map<point, vector> ns);
 
     /*!\brief Construct the domain from filenames.
      * 
@@ -126,7 +126,7 @@
    \mathcal{B}u = g &\text{on } \partial\Omega.
    \end{cases}
    \f]
-   */
+  */
   class BVP{
   public:
     /*! \brief Create a boundary value problem.
@@ -142,29 +142,19 @@
      * \param g_in - An std::map or a realfunc giving the values that the 
      * boundary operator must take.
      */
-    //FIXME: Templatise this later
-    BVP(shared_ptr<const domain> O, 
-	shared_ptr<const diff_op> L_in, 
-	shared_ptr<const bdry_diff_op> B_in, 
-	const realfunc& f_in, 
-	const realfunc& g_in);
-    BVP(shared_ptr<const domain> O, 
-	shared_ptr<const diff_op> L_in, 
-	shared_ptr<const bdry_diff_op> B_in, 
-	const realfunc& f_in, 
-	const map<point,double>& g_in);
+    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 map<point,double>& f_in, 
-	const realfunc& g_in);
-    BVP(shared_ptr<const domain> O, 
-	shared_ptr<const diff_op> L_in, 
-	shared_ptr<const bdry_diff_op> B_in, 
-	const map<point,double>& f_in, 
-	const map<point,double>& g_in);
-
-
+        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);
+    }
 
     virtual ~BVP() {};
 
@@ -215,31 +205,14 @@
   class linear_BVP2 : public BVP{
   public:
     /// Identical to base class constructor.
-    linear_BVP2(shared_ptr<const domain> O, 
-	shared_ptr<const linear_diff_op2> L_in, 
-	shared_ptr<const bdry_diff_op> B_in, 
-	const realfunc& f_in, 
-	const realfunc& g_in)      
-      : BVP(O, L_in, B_in, f_in, g_in){};
-    //FIXME: figure out how to templatise this.
+    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 realfunc& f_in, 
-	const map<point,double>& g_in)
-      : BVP(O, L_in, B_in, f_in, g_in){};
-    linear_BVP2(shared_ptr<const domain> O, 
-	shared_ptr<const linear_diff_op2> L_in, 
-	shared_ptr<const bdry_diff_op> B_in, 
-	const map<point,double>& f_in, 
-	const realfunc& g_in)
-      : BVP(O, L_in, B_in, f_in, g_in){};
-    linear_BVP2(shared_ptr<const domain> O, 
-	shared_ptr<const linear_diff_op2> L_in, 
-	shared_ptr<const bdry_diff_op> B_in, 
-	const map<point,double>& f_in, 
-	const map<point,double>& g_in)
-      : BVP(O, L_in, B_in, f_in, g_in){};
+                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;
--- a/src/include/ddm.hpp
+++ b/src/include/ddm.hpp
@@ -16,7 +16,8 @@
 #include "diff_op.hpp"
 #include "func.hpp"
 
-namespace kwantix{
+namespace kwantix
+{
   using std::map;
   using std::set;
   using boost::shared_ptr;
@@ -30,11 +31,12 @@
    * decomposition, the alternate Bprime operator is applied instead
    * of the B operator.
    */
-  class ddm_bdry_diff_op : public bdry_diff_op{
+  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);
+                     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:
@@ -51,7 +53,7 @@
   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);
+        shared_ptr<const BVP> thebvp);
         
     ///Relative tolerance. Defaults to 1e-5.
     void set_tolerance(double tol);
@@ -84,12 +86,12 @@
   class additive_schwarz_ddm : public ddm{
   public:
     additive_schwarz_ddm(const set<shared_ptr<const domain> >&, 
-			 shared_ptr<const linear_BVP2> bvp);
+                         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;
+        shared_ptr<interpolator<RBF> > > phis;
 
     /*! \brief Say which interpolators are at work at given point.
      *
@@ -104,13 +106,13 @@
 
     ///Evaluate using the average of relevant interpolators.
     double avg_interp(set<shared_ptr<const interpolator<RBF> > > 
-		      relevant_interpolators,
-		      const point& p) const;
+                      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;
     
   };
@@ -135,7 +137,7 @@
    * 
    * 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.
@@ -154,15 +156,15 @@
      * 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);
+                       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);
+                       set<point> bdry, map<point, vector> ns);
       
     /// Give the domains with which this domain overlaps.
     set<shared_ptr<const overlapping_domain> > get_domains() const;
@@ -171,11 +173,11 @@
     shared_ptr<const overlapping_domain> which_domain(const point& p) const; 
     
     friend void set_overlapper_info(set<shared_ptr<overlapping_domain> > 
-				    domains);
+                                    domains);
 
     ///Set overlapper info manually.
     void set_overlapper_info(const point& p, 
-			     const shared_ptr<overlapping_domain> o);
+                             const shared_ptr<overlapping_domain> o);
 
   private:
     /// The domains with which this one overlaps.
@@ -193,7 +195,7 @@
    * with which one.
    */
   void set_overlapper_info(set<shared_ptr<overlapping_domain> > 
-			   domains);
+                           domains);
 }
 
 #endif // __DDM_HPP__
--- a/src/vtkplot.cpp
+++ b/src/vtkplot.cpp
@@ -28,26 +28,23 @@
 #include <string>
 using std::string;
 
-#include "include/linalg.hpp"
+#include "include/vtkplotting.hpp"
 
 typedef double (*Func3d)(double, double);
 
 void InitOffscreen(bool offscreen);
-vtkSmartPointer<vtkDelaunay2D> TriangulateTerrain(Func3d func);
+
+vtkSmartPointer<vtkDelaunay2D> 
+TriangulateTerrain(kwantix::matrix points);
 
 void PlotPoints(vtkSmartPointer<vtkDelaunay2D> delaunay, 
                 bool offscreen=true,
                 const string& filename="vtkscreenshot.png");
 
-double F(double x, double y)
-{
-  return atan(x*y*sin(y));
-}
-
 int run ()
 {
-	auto delaunay = TriangulateTerrain(&F);	
-  PlotPoints(delaunay,true);
+  //	auto delaunay = TriangulateTerrain();	
+  //  PlotPoints(delaunay,true);
 	return 0;
 }