diff --git a/src/mesh.h b/src/mesh.h
--- a/src/mesh.h
+++ b/src/mesh.h
@@ -45,385 +45,384 @@
 using namespace std;
 // new queue which rejects duplicate elements
 template<class T, class C = deque<T> > class unique_queue : public queue<T,C> {
-	
-public:
-	typedef typename C::value_type value_type;
-	// reimplements push: reject element if it exists already
-	void push(const value_type &x) {
-		if (find (queue<T,C>::c.begin(),queue<T,C>::c.end(),x)==queue<T,C>::c.end()) {
-			queue<T,C>::c.push_back(x);
-		}
-	}
-	void clear(void) {
-		queue<T,C>::c.clear();
-	}
+
+ public:
+ typedef typename C::value_type value_type;
+ // reimplements push: reject element if it exists already
+ void push(const value_type &x) {
+   if (find (queue<T,C>::c.begin(),queue<T,C>::c.end(),x)==queue<T,C>::c.end()) {
+     queue<T,C>::c.push_back(x);
+   }
+ }
+ void clear(void) {
+   queue<T,C>::c.clear();
+ }
 };
 
 template<class P> P& deref_ptr ( P *obj) { return *obj; }
 
 
 class Mesh {
-	
-	friend class Cell;
-	friend class Node;
-	friend class FigureEditor;
-	
-public: 
-	Mesh(void) {
-		// Make sure the reserved value is large enough if a cell is added
-		// in "Divide" when the capacity is insufficient, "cells" might be
-		// relocated including the current Cell (i.e. the value of *this)
-		// calling "Mesh::IncreaseCapacityIfNecessary" (from another
-		// object than Cell, e.g. Mesh) before entering Divide will solve
-		// this issue (solved now).
-		cells.reserve(2);
-		nodes.reserve(500);
-		
-		time = 0.;
-		plugin = 0;
-	};
-	~Mesh(void) {
-		delete boundary_polygon;
-	};
-	
-	void Clean(void);
-	Cell &EllipticCell(double xc, double yc, double ra, double rb, int nnodes=10, double rotation=0);
-	Cell &CircularCell(double xc, double yc, double r, int nnodes=10) {
-		return EllipticCell(xc, yc, r, r, nnodes, 0);
-	}
-	Cell &LeafPrimordium(int n, double pet_length);
-	Cell &LeafPrimordium2(int n);
-	Cell *RectangularCell(const Vector ll, const Vector ur, double rotation = 0);
-	void CellFiles(const Vector ll, const Vector ur);
-	
-	inline Cell &getCell(int i) {
-		if ((unsigned)i<cells.size())
-			return *cells[i];
-		else {
-                        #ifdef QDEBUG
-                        qDebug() << i << endl;
-                        qDebug() << "size is " << cells.size() << endl;
-                        #endif
-			abort();
-		}
-	}
-	
-	inline Node &getNode(int i) {
-		return *nodes[i];    
-	}
-	
-	//double Diffusion(void);
-	inline int size(void) {
-		return cells.size();
-	}
-	inline int nnodes(void) {
-		return nodes.size();
-	}
+
+  friend class Cell;
+  friend class Node;
+  friend class FigureEditor;
+
+ public: 
+  Mesh(void) {
+    // Make sure the reserved value is large enough if a cell is added
+    // in "Divide" when the capacity is insufficient, "cells" might be
+    // relocated including the current Cell (i.e. the value of *this)
+    // calling "Mesh::IncreaseCapacityIfNecessary" (from another
+    // object than Cell, e.g. Mesh) before entering Divide will solve
+    // this issue (solved now).
+    cells.reserve(2);
+    nodes.reserve(500);
+
+    time = 0.;
+    plugin = 0;
+  };
+  ~Mesh(void) {
+    delete boundary_polygon;
+  };
+
+  void Clean(void);
+  Cell &EllipticCell(double xc, double yc, double ra, double rb, int nnodes=10, double rotation=0);
+  Cell &CircularCell(double xc, double yc, double r, int nnodes=10) {
+    return EllipticCell(xc, yc, r, r, nnodes, 0);
+  }
+  Cell &LeafPrimordium(int n, double pet_length);
+  Cell &LeafPrimordium2(int n);
+  Cell *RectangularCell(const Vector ll, const Vector ur, double rotation = 0);
+  void CellFiles(const Vector ll, const Vector ur);
+
+  inline Cell &getCell(int i) {
+    if ((unsigned)i<cells.size())
+      return *cells[i];
+    else {
+#ifdef QDEBUG
+      qDebug() << i << endl;
+      qDebug() << "size is " << cells.size() << endl;
+#endif
+      abort();
+    }
+  }
+
+  inline Node &getNode(int i) {
+    return *nodes[i];    
+  }
+
+  //double Diffusion(void);
+  inline int size(void) {
+    return cells.size();
+  }
+  inline int nnodes(void) {
+    return nodes.size();
+  }
+
+  template<class Op> void LoopCells(Op f) {
+    for (vector <Cell *>::iterator i=cells.begin();
+	 i!=cells.end();
+	 i++) {
+      f(**i);
+    }
+  }
+
+  template<class Op> void LoopWalls(Op f) {
+    for (list <Wall *>::iterator i=walls.begin();
+	 i!=walls.end();
+	 i++) {
+      f(**i);
+    }
+  }
 
-	template<class Op> void LoopCells(Op f) {
-		for (vector <Cell *>::iterator i=cells.begin();
-			 i!=cells.end();
-			 i++) {
-			f(**i);
-		}
-	}
-	
-	template<class Op> void LoopWalls(Op f) {
-		for (list <Wall *>::iterator i=walls.begin();
-			 i!=walls.end();
-			 i++) {
-			f(**i);
-		}
-	}
-	
-	// if the amount of cells might increase, during looping, use this template
-	template<class Op> void LoopCurrentCells(Op f) {
-		vector<Cell *> current_cells = cells;
-		for (vector <Cell *>::iterator i=current_cells.begin();
-			 i!=current_cells.end();
-			 i++) {
-			f(**i);
-			
-		}
-	}
-	
-	template<class Op> void LoopNodes(Op f) {
-		for (vector<Node *>::iterator i=nodes.begin();
-			 i!=nodes.end();
-			 i++) {
-			f(**i); 
-		}
-	}
-	
-	template<class Op> void RandomlyLoopNodes(Op f) {
-		
-		MyUrand r(shuffled_nodes.size());
-		random_shuffle(shuffled_nodes.begin(),shuffled_nodes.end(),r);
-		
-		for (vector<Node *>::const_iterator i=shuffled_nodes.begin();
-			 i!=shuffled_nodes.end();
-			 i++) {
-			f(*shuffled_nodes[*i]);
-		}
-		
-	}
-	
-	template<class Op> void RandomlyLoopCells(Op f) {
-		
-		MyUrand r(shuffled_cells.size());
-		random_shuffle(shuffled_cells.begin(),shuffled_cells.end(),r);
-		
-		for (vector<Cell *>::const_iterator i=shuffled_cells.begin();
-			 i!=shuffled_cells.end();
-			 i++) {
-			f(*shuffled_cells[*i]);
-		}
-		
-		
-	}
-	
-	template<class Op1, class Op2> void LoopCells(Op1 f, Op2 &g) {
-		for (vector<Cell *>::iterator i=cells.begin();
-			 i!=cells.end();
-			 i++) {
-			f(**i,g); 
-		}
-	}
-	
-	template<class Op1, class Op2, class Op3> void LoopCells(Op1 f, Op2 &g, Op3 &h) {
-		for (vector<Cell *>::iterator i=cells.begin();
-			 i!=cells.end();
-			 i++) {
-			f(**i,g,h); 
-		}
-	}
-	
-	void DoCellHouseKeeping(void) {
-		vector<Cell *> current_cells = cells;
-		for (vector<Cell *>::iterator i = current_cells.begin();
-			 i != current_cells.end();
-			 i ++) {
-			plugin->CellHouseKeeping(*i);
-			
-			// Call functions of Cell that cannot be called from CellBase, including Division
-			if ((*i)->flag_for_divide) {
-				if ((*i)->division_axis) {
-					(*i)->DivideOverAxis(*(*i)->division_axis);
-					delete (*i)->division_axis;
-					(*i)->division_axis = 0;
-				} else {
-					(*i)->Divide();
-				}
-				(*i)->flag_for_divide=false;
-			}
-		}
-	}
-	
-	// Apply "f" to cell i
-	// i.e. this is an adapter which allows you to call a function
-	// operating on Cell on its numeric index index
-	template<class Op> void cell_index_adapter(Op f,int i) {
-		f(cells[i]);
-	}
-	
-	double DisplaceNodes(void);
-	
-	void BoundingBox(Vector &LowerLeft, Vector &UpperRight);
-	int NEqs(void) {     int nwalls = walls.size();
-		int ncells =cells.size();
-		int nchems = Cell::NChem();
-		
-		// two eqs per chemical for each walls, and one eq per chemical for each cell
-		// This is for generality. For a specific model you may optimize
-		// this by removing superfluous (empty) equations.
-		int neqs = 2 * nwalls * nchems + ncells * nchems;
-		
-		return neqs;
+  // if the amount of cells might increase, during looping, use this template
+  template<class Op> void LoopCurrentCells(Op f) {
+    vector<Cell *> current_cells = cells;
+    for (vector <Cell *>::iterator i=current_cells.begin();
+	 i!=current_cells.end();
+	 i++) {
+      f(**i);
+
+    }
+  }
+
+  template<class Op> void LoopNodes(Op f) {
+    for (vector<Node *>::iterator i=nodes.begin();
+	 i!=nodes.end();
+	 i++) {
+      f(**i); 
+    }
+  }
+
+  template<class Op> void RandomlyLoopNodes(Op f) {
+
+    MyUrand r(shuffled_nodes.size());
+    random_shuffle(shuffled_nodes.begin(),shuffled_nodes.end(),r);
+
+    for (vector<Node *>::const_iterator i=shuffled_nodes.begin();
+	 i!=shuffled_nodes.end();
+	 i++) {
+      f(*shuffled_nodes[*i]);
+    }
+  }
+
+  template<class Op> void RandomlyLoopCells(Op f) {
+
+    MyUrand r(shuffled_cells.size());
+    random_shuffle(shuffled_cells.begin(),shuffled_cells.end(),r);
+
+    for (vector<Cell *>::const_iterator i=shuffled_cells.begin();
+	 i!=shuffled_cells.end();
+	 i++) {
+      f(*shuffled_cells[*i]);
+    }
+  }
+
+  template<class Op1, class Op2> void LoopCells(Op1 f, Op2 &g) {
+    for (vector<Cell *>::iterator i=cells.begin();
+	 i!=cells.end();
+	 i++) {
+      f(**i,g); 
+    }
+  }
+
+  template<class Op1, class Op2, class Op3> void LoopCells(Op1 f, Op2 &g, Op3 &h) {
+    for (vector<Cell *>::iterator i=cells.begin();
+	 i!=cells.end();
+	 i++) {
+      f(**i,g,h); 
+    }
+  }
+
+  void DoCellHouseKeeping(void) {
+    vector<Cell *> current_cells = cells;
+    for (vector<Cell *>::iterator i = current_cells.begin();
+	 i != current_cells.end();
+	 i ++) {
+      plugin->CellHouseKeeping(*i);
+
+      // Call functions of Cell that cannot be called from CellBase, including Division
+      if ((*i)->flag_for_divide) {
+	if ((*i)->division_axis) {
+	  (*i)->DivideOverAxis(*(*i)->division_axis);
+	  delete (*i)->division_axis;
+	  (*i)->division_axis = 0;
+	} else {
+	  (*i)->Divide();
 	}
-	void IncreaseCellCapacityIfNecessary(void) {
-		
-		return;
-		// cerr << "Entering Mesh::IncreaseCellCapacityIfNecessary \n";
-		// make sure we always have enough space 
-		// to have each cell divide at least once
-		//
-		// Note that we must do this, because Cell::Divide pushes a new Cell
-		// onto Mesh::cells. As a result, Mesh::cells might be relocated 
-		// if we are _within_ a Cell object: i.e. pointer "this" will be changed!!
-		// 
-		// An alternative solution could be to make "Mesh::cells" a list,
-		// but this won't work because we need random access for 
-		// the Monte Carlo algorithm.
-		
-		if (2*cells.size()>cells.capacity()) {
-			cerr << "Increasing capacity to "  << 2*cells.capacity() << endl;
-			cerr << "Current capacity is " << cells.capacity() << endl;
-			cells.reserve(cells.capacity()*2);
-		}
-	}
-	
-	void ReserveMoreCells(int n) {
-		if (nodes.size()+n>nodes.capacity()) {
-			nodes.reserve(size()+n);
-		}
-	}
-	double Area(void);
-	double MeanArea(void) {
-		double sum=0.;
-		for (vector<Cell *>::const_iterator i=cells.begin();
-			 i!=cells.end();
-			 i++) {
-			sum+=(*i)->Area();
-		}
-		return sum/(double)NCells();
-	}
-	
-	void SetBaseArea(void);
-	int NCells(void) const {
-		return cells.size();
-	}
-	inline int NNodes(void) const {
-		return nodes.size();
-	}
-	void PrintQueue(ostream &os) {
-		while (!node_insertion_queue.empty()) {
-			os << node_insertion_queue.front() << endl;
-			node_insertion_queue.pop();
-		}
-	}
-	
-	void InsertNodes(void) {
-		// insert the nodes in the insertion queue
-		while (!node_insertion_queue.empty()) {
-			
-			//cerr << node_insertion_queue.front() << endl;
-			InsertNode(node_insertion_queue.front());
-			node_insertion_queue.pop();
-		}
-		
-	}
-	
-	void Clear(); 
-	
-	void ReactDiffuse( double delta_t = 1 );
-	double SumChemical(int ch);
-	void SetChemical(int ch, double value) {
-		for (vector<Cell *>::iterator c=cells.begin();
-			 c!=cells.end();
-			 c++) {
-			(*c)->chem[ch]=value;
-		}
-	}
-	
-	// used for interacing with ODE-solvers (e.g. NRCRungeKutta)
-	void setValues(double x, double *y);
-	double *getValues(int *neqs);
-	void Derivatives(double *derivs);
+	(*i)->flag_for_divide=false;
+      }
+    }
+  }
+
+  // Apply "f" to cell i
+  // i.e. this is an adapter which allows you to call a function
+  // operating on Cell on its numeric index index
+  template<class Op> void cell_index_adapter(Op f,int i) {
+    f(cells[i]);
+  }
+
+  double DisplaceNodes(void);
+
+  void BoundingBox(Vector &LowerLeft, Vector &UpperRight);
+  int NEqs(void) {     int nwalls = walls.size();
+    int ncells =cells.size();
+    int nchems = Cell::NChem();
+
+    // two eqs per chemical for each walls, and one eq per chemical for each cell
+    // This is for generality. For a specific model you may optimize
+    // this by removing superfluous (empty) equations.
+    int neqs = 2 * nwalls * nchems + ncells * nchems;
+
+    return neqs;
+  }
+  void IncreaseCellCapacityIfNecessary(void) {
+
+    return;
+    // cerr << "Entering Mesh::IncreaseCellCapacityIfNecessary \n";
+    // make sure we always have enough space 
+    // to have each cell divide at least once
+    //
+    // Note that we must do this, because Cell::Divide pushes a new Cell
+    // onto Mesh::cells. As a result, Mesh::cells might be relocated 
+    // if we are _within_ a Cell object: i.e. pointer "this" will be changed!!
+    // 
+    // An alternative solution could be to make "Mesh::cells" a list,
+    // but this won't work because we need random access for 
+    // the Monte Carlo algorithm.
+
+    if (2*cells.size()>cells.capacity()) {
+      cerr << "Increasing capacity to "  << 2*cells.capacity() << endl;
+      cerr << "Current capacity is " << cells.capacity() << endl;
+      cells.reserve(cells.capacity()*2);
+    }
+  }
+
+  void ReserveMoreCells(int n) {
+    if (nodes.size()+n>nodes.capacity()) {
+      nodes.reserve(size()+n);
+    }
+  }
+  double Area(void);
+  double MeanArea(void) {
+    double sum=0.;
+    for (vector<Cell *>::const_iterator i=cells.begin();
+	 i!=cells.end();
+	 i++) {
+      sum+=(*i)->Area();
+    }
+    return sum/(double)NCells();
+  }
+
+  void SetBaseArea(void);
+  int NCells(void) const {
+    return cells.size();
+  }
+  inline int NNodes(void) const {
+    return nodes.size();
+  }
+  void PrintQueue(ostream &os) {
+    while (!node_insertion_queue.empty()) {
+      os << node_insertion_queue.front() << endl;
+      node_insertion_queue.pop();
+    }
+  }
+
+  void InsertNodes(void) {
+    // insert the nodes in the insertion queue
+    while (!node_insertion_queue.empty()) {
+
+      //cerr << node_insertion_queue.front() << endl;
+      InsertNode(node_insertion_queue.front());
+      node_insertion_queue.pop();
+    }
+
+  }
+
+  void Clear(); 
+
+  void ReactDiffuse( double delta_t = 1 );
+  double SumChemical(int ch);
+  void SetChemical(int ch, double value) {
+    for (vector<Cell *>::iterator c=cells.begin();
+	 c!=cells.end();
+	 c++) {
+      (*c)->chem[ch]=value;
+    }
+  }
+
+  // used for interacing with ODE-solvers (e.g. NRCRungeKutta)
+  void setValues(double x, double *y);
+  double *getValues(int *neqs);
+  void Derivatives(double *derivs);
 #ifdef QTGRAPHICS
-	inline void DrawBoundary(QGraphicsScene *c) {
-		boundary_polygon->Draw(c);
-	}
-	void DrawNodes(QGraphicsScene *c) const;
-	
+  inline void DrawBoundary(QGraphicsScene *c) {
+    boundary_polygon->Draw(c);
+  }
+  void DrawNodes(QGraphicsScene *c) const;
+
 #endif
-	double max_chem;
-	
-	void XMLSave(const char *docname, xmlNode *settings=0) const;
-	void XMLRead(const char *docname, xmlNode **settings=0, bool geometry = true, bool pars = true, bool simtime = true);
-	void XMLReadPars(const xmlNode * root_node);
-	void XMLReadGeometry(const xmlNode *root_node);
-	void XMLReadSimtime(const xmlNode *root_node);
-	void XMLReadNodes(xmlNode *cur);
-	void XMLReadCells(xmlNode *cur);
-	void XMLParseTree(const xmlNode * root_node);
-	void XMLReadWalls(xmlNode *cur, vector<Wall *> *tmp_cells);
-	void XMLReadWallsToCells(xmlNode *root, vector<Wall *> *tmp_walls);
-	void XMLReadNodeSets(xmlNode *root);
-	void XMLReadNodeSetsToNodes(xmlNode *root);
-	void PerturbChem(int chemnum, double range);
-	void CleanUpCellNodeLists(void);
-	void CleanUpWalls(void);
-	void CutAwayBelowLine( Vector startpoint, Vector endpoint );
-	void CutAwaySAM(void);
-	void RepairBoundaryPolygon(void);
-	void Rotate(double angle, Vector center);
-	void PrintWallList( void );
-	void TestIllegalWalls(void);
-	Vector FirstConcMoment(int chem);
-	inline Vector Centroid(void) {
-		return boundary_polygon->Centroid();
-	}
-	
-	inline Vector Offset(void) {
-		return boundary_polygon->Offset();
-	}
-	
-	inline double Factor(void) {
-		return boundary_polygon->Factor();
-	}
-	
-	void DeleteLooseWalls(void);
-	void FitLeafToCanvas(double width, double height);
-	void AddNodeSet(NodeSet *node_set) {
-		node_sets.push_back(node_set);
-	}
-	
-	void CleanChemicals(const vector<double> &clean_chem);
-	void CleanTransporters(const vector<double> &clean_transporters);
-	void RandomizeChemicals(const vector<double> &max_chem, const vector<double> &max_transporters);
-	inline double getTime(void) const { return time; }
-	string getTimeHours(void) const; 
-	inline void setTime(double t) { time = t; }
-	double CalcProtCellsWalls(int ch) const;  
-	void SettoInitVals(void);
-	QVector<qreal> VertexAngles(void);
-	QVector< QPair<qreal,int> > VertexAnglesValues(void);
-	void SetSimPlugin(SimPluginInterface *new_plugin) {
-		plugin=new_plugin;
-	}
-	QString ModelID(void) { return plugin?plugin->ModelID():QString("undefined"); }
-	void StandardInit(void);	
+  double max_chem;
+
+  void XMLSave(const char *docname, xmlNode *settings=0) const;
+  void XMLRead(const char *docname, xmlNode **settings=0, bool geometry = true, bool pars = true, bool simtime = true);
+  void XMLReadPars(const xmlNode * root_node);
+  void XMLReadGeometry(const xmlNode *root_node);
+  void XMLReadSimtime(const xmlNode *root_node);
+  void XMLReadNodes(xmlNode *cur);
+  void XMLReadCells(xmlNode *cur);
+  void XMLParseTree(const xmlNode * root_node);
+  void XMLReadWalls(xmlNode *cur, vector<Wall *> *tmp_cells);
+  void XMLReadWallsToCells(xmlNode *root, vector<Wall *> *tmp_walls);
+  void XMLReadNodeSets(xmlNode *root);
+  void XMLReadNodeSetsToNodes(xmlNode *root);
+  void PerturbChem(int chemnum, double range);
+  void CleanUpCellNodeLists(void);
+  void CleanUpWalls(void);
+  void CutAwayBelowLine( Vector startpoint, Vector endpoint );
+  void CutAwaySAM(void);
+  void RepairBoundaryPolygon(void);
+  void Rotate(double angle, Vector center);
+  void PrintWallList( void );
+  void TestIllegalWalls(void);
+  Vector FirstConcMoment(int chem);
+  inline Vector Centroid(void) {
+    return boundary_polygon->Centroid();
+  }
+
+  inline Vector Offset(void) {
+    return boundary_polygon->Offset();
+  }
+
+  inline double Factor(void) {
+    return boundary_polygon->Factor();
+  }
+
+  void DeleteLooseWalls(void);
+  void FitLeafToCanvas(double width, double height);
+  void AddNodeSet(NodeSet *node_set) {
+    node_sets.push_back(node_set);
+  }
 
-	Node* findNextBoundaryNode(Node*);
+  void CleanChemicals(const vector<double> &clean_chem);
+  void CleanTransporters(const vector<double> &clean_transporters);
+  void RandomizeChemicals(const vector<double> &max_chem, const vector<double> &max_transporters);
+  inline double getTime(void) const { return time; }
+  string getTimeHours(void) const; 
+  inline void setTime(double t) { time = t; }
+  double CalcProtCellsWalls(int ch) const;  
+  void SettoInitVals(void);
+  QVector<qreal> VertexAngles(void);
+  QVector< QPair<qreal,int> > VertexAnglesValues(void);
+  void SetSimPlugin(SimPluginInterface *new_plugin) {
+    plugin=new_plugin;
+  }
+  QString ModelID(void) { return plugin?plugin->ModelID():QString("undefined"); }
+  void StandardInit(void);	
+
+  Node* findNextBoundaryNode(Node*);
+
+ private:
 
-private:
-	
-	// Data members
-	vector<Cell *> cells;
-	vector<Node *> nodes;
-	list<Wall *> walls; // we need to erase elements from this container frequently, hence a list.
-public:
-	vector<NodeSet *> node_sets;
-private:
-	vector<Node *> shuffled_nodes;
-	vector<Cell *> shuffled_cells;
-	unique_queue<Edge> node_insertion_queue;
-	BoundaryPolygon *boundary_polygon;
-	double time;
-	SimPluginInterface *plugin;
-	
-	// Private member functions
-	void AddNodeToCell(Cell *c, Node *n, Node *nb1 , Node *nb2);
-	void AddNodeToCellAtIndex(Cell *c, Node *n, Node *nb1 , Node *nb2, list<Node *>::iterator ins_pos);
-	void InsertNode(Edge &e);
-	inline Node *AddNode(Node *n) {
-		nodes.push_back(n);
-		shuffled_nodes.push_back(n);
-		n->m=this;
-		return n;
-	}
-	
-	inline Cell *AddCell(Cell *c) {
-		cells.push_back(c);
-		shuffled_cells.push_back(c);
-		//cerr << "Shuffled cell indices:  ";
-		/*copy(shuffled_cells.begin(),shuffled_cells.end(),ostream_iterator<int>(cerr," "));
-		 cerr << endl;*/
-		c->m=this;
-		return c;
-	}
-	
-	void CircumCircle(double x1,double y1,double x2,double y2,double x3,double y3,
-					  double *xc,double *yc,double *r);
+  // Data members
+  vector<Cell *> cells;
+  vector<Node *> nodes;
+  list<Wall *> walls; // we need to erase elements from this container frequently, hence a list.
+ public:
+  vector<NodeSet *> node_sets;
+ private:
+  vector<Node *> shuffled_nodes;
+  vector<Cell *> shuffled_cells;
+  unique_queue<Edge> node_insertion_queue;
+  BoundaryPolygon *boundary_polygon;
+  double time;
+  SimPluginInterface *plugin;
+
+  // Private member functions
+  void AddNodeToCell(Cell *c, Node *n, Node *nb1 , Node *nb2);
+  void AddNodeToCellAtIndex(Cell *c, Node *n, Node *nb1 , Node *nb2, list<Node *>::iterator ins_pos);
+  void InsertNode(Edge &e);
+  inline Node *AddNode(Node *n) {
+    nodes.push_back(n);
+    shuffled_nodes.push_back(n);
+    n->m=this;
+    return n;
+  }
+
+  inline Cell *AddCell(Cell *c) {
+    cells.push_back(c);
+    shuffled_cells.push_back(c);
+    //cerr << "Shuffled cell indices:  ";
+    /*copy(shuffled_cells.begin(),shuffled_cells.end(),ostream_iterator<int>(cerr," "));
+      cerr << endl;*/
+    c->m=this;
+    return c;
+  }
+
+  void CircumCircle(double x1,double y1,double x2,double y2,double x3,double y3,
+		    double *xc,double *yc,double *r);
 };
 #endif
+
+/* finis */