EI/VirtualLeaf Changeset - d9ce329d84ae · Centrum Wiskunde & Informatica (CWI)

Changeset - d9ce329d84ae

Parent rev.

Child rev.

[Not reviewed]

default

0 13 0

Michael Guravage - 15 years ago 2010-06-18 09:37:17
michael.guravage@cwi.nl

Made additional formatting changes.
--
user: Michael Guravage <michael.guravage@cwi.nl>
branch 'default'
changed src/VirtualLeaf.cpp
changed src/canvas.cpp
changed src/cell.cpp
changed src/cellbase.cpp
changed src/mesh.cpp
changed src/modelcatalogue.h
changed src/modelelement.h
changed src/node.cpp
changed src/parameter.cpp
changed src/parameter.h
changed src/pardialog.cpp
changed src/pardialog.h
changed src/xmlwrite.cpp

13 files changed with 172 insertions and 543 deletions:

src/VirtualLeaf.cpp

src/canvas.cpp

src/cell.cpp

154

src/cellbase.cpp

src/mesh.cpp

309

src/modelcatalogue.h

src/modelelement.h

src/node.cpp

src/parameter.cpp

src/parameter.h

src/pardialog.cpp

src/pardialog.h

src/xmlwrite.cpp

0 comments (0 inline, 0 general)

src/VirtualLeaf.cpp

➞

Show inline comments

@@ @@ -4,195 +4,193 @@ @@
+ *
  *  The Virtual Leaf is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
+ *
  *  The Virtual Leaf is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU General Public License for more details.
+ *
  *  You should have received a copy of the GNU General Public License
  *  along with the Virtual Leaf.  If not, see <http://www.gnu.org/licenses/>.
+ *
  *  Copyright 2010 Roeland Merks.
+ *
  */
 #include <string>
 #include <fstream>
 #include <sstream>
 #include <cstring>
 #include <functional>
 #include <getopt.h>
 #include <cerrno>
 #include "mesh.h"
 #include "parameter.h"
 #include "random.h"
 #include "pi.h"
 #include "cellitem.h"
 #include "canvas.h"
 #include "cell.h"
 #include "output.h"
 #include <qwidget.h>
 #include <q3process.h>
 #include <qapplication.h>
 #include <QDesktopWidget>
 #include <QGraphicsScene>
 #include <QMessageBox>
 //Added by qt3to4:
 #include <QMouseEvent>
 #include <unistd.h>
 #include <q3textstream.h>
 #ifdef HAVE_QWT
 #include "data_plot.h"
 #endif
 #include <QPalette>
 #include <QBrush>
 #include <QToolTip>
 #include "simplugin.h"
 #include <QPluginLoader>
 #include <QDir>
 #include "modelcatalogue.h"
 static const std::string _module_id("$Id$");
 extern Parameter par;
 MainBase *main_window = 0;
 #ifdef XFIGGRAPHICS
 #define TIMESTEP double Graphics::TimeStep(void)
 #endif
 class PrintNode {
 public:
   void operator() (const Node &n) const
+  {
     cerr << n.Index() << ": " << n <<  endl;
+  }
 };
 class EdgeSource {
 public:
   void operator() (Cell &c) {
     if (c.AtBoundaryP()) {
       cerr << "Cell " << c.Index() << " is a source cell.\n";
       c.SetSource(0,par.source);
     } else {
       cerr << "Cell " << c.Index() << " is _not_ a source cell.\n";
+    }
+  }
 };
 class CellInfo {
 public:
   void operator() (Cell &c,std::ostream &os) const {
     os << "Cell " << c.index << " says: " << endl;
     os << "c.nodes.size() = " << c.nodes.size() << endl;
     for (list<Node *>::iterator i=c.nodes.begin();
 	 i!=c.nodes.end();
 	 i++) {
     for (list<Node *>::iterator i=c.nodes.begin(); i!=c.nodes.end(); i++) {
       cerr << (*i)->Index() << " ";
+    }
     cerr << endl;
+  }
 };
 double PINSum(Cell &c) {
   return c.Chemical(1) + c.SumTransporters(1);// + c.ReduceCellAndWalls<double>( complex_PijAj );
+}
 class DrawCell {
 public:
   void operator() (Cell &c,QGraphicsScene &canvas, MainBase &m) const {
     if (m.ShowBorderCellsP() || c.Boundary()==Cell::None) {
       if (!m.ShowBoundaryOnlyP() && !m.HideCellsP()) {
 	if (m.ShowToolTipsP()) {
 	  QString info_string=QString("Cell %1, chemicals: ( %2, %3, %4, %5, %6)\n %7 of PIN1 at walls.\n Area is %8\n PIN sum is %9\n Circumference is %10\n Boundary type is %11").arg(c.Index()).arg(c.Chemical(0)).arg(c.Chemical(1)).arg(c.Chemical(2)).arg(c.Chemical(3)).arg(c.Chemical(4)).arg(c.SumTransporters(1)).arg(c.Area()).arg(PINSum(c)).arg(c.Circumference()).arg(c.BoundaryStr());
 	  info_string += "\n" + c.printednodelist();
 	  c.Draw(&canvas, info_string);
 	} else {
 	  c.Draw(&canvas);
+	}
+      }
       if (m.ShowCentersP()){
 	c.DrawCenter(&canvas);
+      }
       if (m.ShowFluxesP()){
 	c.DrawFluxes(&canvas, par.arrowsize);
+      }
+    }
+  }
 };
 Mesh mesh;
 bool batch=false;
 void MainBase::Plot(int resize_stride)
+{
   clear();
   static int count=0;
   if (resize_stride) {
     if ( !((++count)%resize_stride) ) {
       FitLeafToCanvas();
+    }
+  }
   mesh.LoopCells(DrawCell(),canvas,*this);
   if (ShowNodeNumbersP())
     mesh.LoopNodes( bind2nd (mem_fun_ref ( &Node::DrawIndex), &canvas ) ) ;
   if (ShowCellNumbersP())
     mesh.LoopCells( bind2nd (mem_fun_ref ( &Cell::DrawIndex), &canvas ) ) ;
   if (ShowCellAxesP())
     mesh.LoopCells( bind2nd (mem_fun_ref ( &Cell::DrawAxis), &canvas ) );
   if (ShowCellStrainP())
     mesh.LoopCells( bind2nd (mem_fun_ref ( &Cell::DrawStrain), &canvas ) );
   if (ShowWallsP())
     mesh.LoopWalls( bind2nd( mem_fun_ref( &Wall::Draw ), &canvas ) );
   if (ShowApoplastsP())
     mesh.LoopWalls( bind2nd( mem_fun_ref( &Wall::DrawApoplast ), &canvas ) );
   if (ShowMeshP())
     mesh.DrawNodes(&canvas);
   if (ShowBoundaryOnlyP())
     mesh.DrawBoundary(&canvas);
   if ( ( batch || MovieFramesP() )) {
     static int frame = 0;
     // frame numbers are sequential for the most frequently written file type.
     // for the less frequently written file type they match the other type
     if (!(count%par.storage_stride) )  {
       stringstream fname;
       fname << par.datadir << "/leaf.";
       fname.fill('0');
       fname.width(6);
       fname << frame << ".jpg";
       if (par.storage_stride <= par.xml_storage_stride) {
 	frame++;
+      }
       // Write high-res JPG snapshot every plot step
       Save(fname.str().c_str(), "JPEG",1024,768);
+    }
     if (!(count%par.xml_storage_stride)) {

src/canvas.cpp

➞

Show inline comments

@@ @@ -202,244 +202,236 @@ void FigureEditor::mousePressEvent(QMous @@
 	  (dynamic_cast<WallItem *>(*it))->OnClick(e);
+	}
+      }
+  }
   FullRedraw();
   moving = 0;
+}
 void FigureEditor::mouseMoveEvent(QMouseEvent* e)
+{
   // User choose "rotation mode" and we can rotate the object around its center of mass
   if (dynamic_cast<Main *>(parent())->RotationModeP()) {
     QPointF p = mapToScene(e->pos());
     p.setX(p.x() / Cell::Magnification());
     p.setY(p.y() / Cell::Magnification());
     // get object's center of mass
     QPointF rotation_midpoint = mesh.Centroid()*Cell::Factor() - Cell::Offset();
     // calculate rotation angle
     double dy = (rotation_midpoint.y() - p.y());
     double dx = (rotation_midpoint.x() - p.x());
     double new_rot_angle = atan2(dx, dy);
     double d_alpha = new_rot_angle - rot_angle;
     rot_angle = new_rot_angle;
     mesh.Rotate(d_alpha, ( Vector(rotation_midpoint) + Cell::Offset() ) / Cell::Factor() );
     dynamic_cast<Main *>(parent())->Plot(0);
     FullRedraw();
     return;
+  }
   if ( moving ) {
     QPointF p = mapToScene(e->pos());
     moving->userMove(p.x() - moving_start.x(),
 		     p.y() - moving_start.y());
     moving_start = p;
     scene()->update();
+  }
   //cerr << "event";
   // keep track of intersection line to interactively cut a growing leaf
   if ( intersection_line ) {
     QPointF sp = intersection_line -> line().p1(); // startpoint
     QPointF ep = mapToScene(e->pos()); // endpoint
     intersection_line -> setLine( QLineF(sp, ep) );
     scene()->update();
     // Need this for Mac
     FullRedraw();
+  }
+}
 //void FigureEditor::contentsMouseReleaseEvent(QMouseEvent* e)
 void FigureEditor::mouseReleaseEvent(QMouseEvent* e)
+{
   emit MouseReleased();
   // intersection line for leaf was finished now.
   if (e->button()==Qt::LeftButton) {
     if (intersection_line ) {
 #ifdef QDEBUG
       qDebug() << "Trying to cut leaf" << endl;
 #endif
       QPointF sp = intersection_line -> line().p1(); // startpoint
       QPointF ep = mapToScene(e->pos());
       intersection_line -> setLine( QLineF(sp, ep) );
       intersection_line -> show();
       vector <CellItem *> intersected_cells = getIntersectedCells();
       // no cells selected, do nothing
       if (intersected_cells.size()==0) {
 #ifdef QDEBUG
 	qDebug() << "No cells detected :-(" << endl;
 #endif
 	return;
+      }
       Vector startpoint = Vector(sp.x(), sp.y()) / Cell::Factor() - Cell::Offset();
       Vector endpoint = Vector(ep.x(), ep.y()) / Cell::Factor() - Cell::Offset();
       NodeSet *node_set = new NodeSet;
       for (vector<CellItem *>::iterator it = intersected_cells.begin();
 	   it != intersected_cells.end();
 	   it++) {
       for (vector<CellItem *>::iterator it = intersected_cells.begin(); it != intersected_cells.end(); it++) {
 	(*it)->setBrush(QBrush("purple"));
 	Cell &c=(*it)->getCell();
 	// sometimes the cell hasn't properly divided yet before the
 	// next division is called?  so check for it?  let's find a way
 	// to do this later. Note that this functionality currently
 	// might result in a segmentation fault for users who are
 	// quickly dragging and releasing division lines...
 	scene()->update();
 #ifdef QDEBUG
 	qDebug() << "Dividing Cell " << c.Index() << endl;
 #endif
 	c.DivideOverGivenLine( startpoint, endpoint, true, node_set);
+      }
       node_set->CleanUp();
       mesh.AddNodeSet(node_set);
 #ifdef QDEBUG
       qDebug() << "Done DivideOverGivenLine" << endl;
 #endif
       mesh.TestIllegalWalls();
       // Do the actual cutting and removing
       if (intersected_cells.size()) {
 	mesh.CutAwayBelowLine( startpoint, endpoint );
 	// Correct flags of nodeset
 	for (
 	     NodeSet::iterator i = node_set->begin();
 	     i != node_set->end();
 	     i++) {
 	for (NodeSet::iterator i = node_set->begin(); i != node_set->end(); i++) {
 	  (*i)->SetSAM();
 	  (*i)->SetBoundary();
+	}
 	// Make cells attached to nodeset part of the boundary
 	// This is a temporary solution for the following:
 	//   If a cell attached to a NodeSet divides (with a division plane intersecting the node set),
 	//   we must insert a new node into the node set.
 	// For now, we add a layer of "virtual cells" inbetween.
 	list<Cell *> cells_attached_to_nodeset = node_set->getCells();
 	for ( list<Cell *>::iterator c = cells_attached_to_nodeset.begin();
 	      c != cells_attached_to_nodeset.end();
 	      c++) {
 	for ( list<Cell *>::iterator c = cells_attached_to_nodeset.begin(); c != cells_attached_to_nodeset.end(); c++) {
 	  (*c)->SetBoundary(Cell::SAM);
+	}
 #ifdef QDEBUG
 	qDebug() << "Done CutAwayBelowLine" << endl;
 #endif
 	mesh.TestIllegalWalls();
 	mesh.RepairBoundaryPolygon();
 #ifdef QDEBUG
 	qDebug() << "Done RepairBoundaryPolygon" << endl;
 #endif
 	mesh.TestIllegalWalls();
 	mesh.CleanUpWalls();
 #ifdef QDEBUG
 	qDebug() << "Done CleanUpWalls" << endl;
 #endif
 	mesh.TestIllegalWalls();
+      }
       dynamic_cast<Main *>(parent())->Plot();
 #ifdef QDEBUG
       qDebug() << "NodeSet of cutting line: " << *node_set << endl;
 #endif
+    }
   } else
     if (e->button()==Qt::RightButton) {
       if (intersection_line /* && !angle_line*/) {
 	QPointF p = mapToScene(e->pos());
 	QPointF sp = intersection_line->line().p1();
 	viewport()->setMouseTracking( TRUE );
+      }
+    }
+}
 // returns a vector of pointer to cells colliding with intersection line
 vector <CellItem *> FigureEditor::getIntersectedCells(void)
+{
   vector <CellItem *> colliding_cells;
   QList<QGraphicsItem *> l = intersection_line->collidingItems( );
 #ifdef QDEBUG
   qDebug() <<  "l.size() = " << l.size() << endl;
 #endif
   for (QList<QGraphicsItem *>::Iterator it=l.begin(); it!=l.end(); ++it) {
 #ifdef QDEBUG
     qDebug() << typeid(**it).name() << endl;
 #endif
     if ( !strcmp(typeid(**it).name(),"8CellItem") ) {
       colliding_cells.push_back(dynamic_cast<CellItem *>(*it));
+    }
+  }
   delete intersection_line;
   intersection_line = 0;
   return colliding_cells;
+}
 void FigureEditor::FullRedraw(void)
+{
   QList<QRectF> rl;
   rl << sceneRect();
   updateScene(rl);
+}
 NodeItem *FigureEditor::selectedNodeItem(QList<QGraphicsItem *> graphicItems) const
+{
   NodeItem *nodeItem = 0;
   // graphicItems is a list of the QgraphicItems under the mouse click event
   QList<QGraphicsItem *>::Iterator it = graphicItems.begin();
   for (; it != graphicItems.end(); ++it) {
     if ( !strcmp(typeid(**it).name(),"8NodeItem")) {
       // the object under the mouse click is a Nodeitem
       nodeItem = dynamic_cast<NodeItem *>(*it);
       // indicate we've selected it
       nodeItem->setBrush(dark_red);
       break;
+    }
+  }
   return nodeItem;
+}
 void FigureEditor::insertNode(QPointF p)
+{

src/cell.cpp

➞

Show inline comments

 /*
+ *
  *  This file is part of the Virtual Leaf.
+ *
  *  The Virtual Leaf is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
+ *
  *  The Virtual Leaf is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU General Public License for more details.
+ *
  *  You should have received a copy of the GNU General Public License
  *  along with the Virtual Leaf.  If not, see <http://www.gnu.org/licenses/>.
+ *
  *  Copyright 2010 Roeland Merks.
+ *
  */
 #include <QDebug>
 #include <string>
 #include "cell.h"
 #include "node.h"
 #include "mesh.h"
 #include "tiny.h"
 #include "nodeset.h"
 #include "cellitem.h"
 #include "nodeitem.h"
 #include "qcanvasarrow.h"
 #include "parameter.h"
 static const std::string _module_id("$Id$");
 extern Parameter par;
 double Cell::factor=1.;
 double Cell::offset[3]={0,0,0};
 Cell::Cell(void) : CellBase()
+{
   m=0;
+}
 Cell::Cell(double x, double y, double z) : CellBase(x,y,z)
+{
   m=0;
+}
 Cell::Cell(const Cell &src) :  CellBase(src)
+{
   m=src.m;
+}
 bool Cell::Cmp(Cell *c) const { return this->Index() < c->Index(); }
 bool Cell::Eq(Cell *c) const { return this->Index() == c->Index(); }
 Cell Cell::operator=(const Cell &src)
+{
   CellBase::operator=(src);
   m=src.m;
   return *this;
+}
 //Cell(void) : CellBase() {}
 void Cell::DivideOverAxis(Vector axis)
+{
   // Build a wall
   // ->  find the position of the wall
   // better look for intersection with a simple line intersection algorithm as below?
   // this leads to some exceptions: e.g. dividing a horizontal rectangle.
   // leaving it like this for the time being
   if (dead) return;
   Vector centroid=Centroid();
   double prev_cross_z=(axis * (centroid - *(nodes.back()) ) ).z ;
   ItList new_node_locations;
   for (list<Node *>::iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (list<Node *>::iterator i=nodes.begin(); i!=nodes.end(); i++) {
     // cross product to detect position of division
     Vector cross = axis * (centroid - *(*i));
     if (cross.z * prev_cross_z < 0 ) {
       new_node_locations.push_back(i);
+    }
     prev_cross_z=cross.z;
+  }
   DivideWalls(new_node_locations, centroid, centroid+axis);
+}
 double Cell::MeanArea(void)
+{
   return m->MeanArea();
+}
 void Cell::Apoptose(void)
+{
   // First kill walls
 #ifdef QDEBUG
   qDebug() << "This is cell " << Index() << endl;
   qDebug() << "Number of walls: " << walls.size() << endl;
 #endif
   for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
 #ifdef QDEBUG
     qDebug() << "Before apoptosis, wall " << (*w)->Index() << " says: c1 = "
 	     << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl;
 #endif
+  }
   for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
   for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
     bool illegal_flag = false;
     if ((*w)->c1 == (*w)->c2 )
       illegal_flag=true;
     if ((*w)->c1 == this) {
       // invert wall?
       (*w)->c1 = (*w)->c2;
       (*w)->c2 = m->boundary_polygon;
       Node *n1 = (*w)->n1;
       (*w)->n1 = (*w)->n2;
       (*w)->n2 = n1;
     } else {
       (*w)->c2 = m->boundary_polygon;
+    }
 #ifdef QDEBUG
     if (illegal_flag && (*w)->c1==(*w)->c2) {
       qDebug() << "I created an illegal wall." << endl;
+    }
 #endif
     if ( ((*w)->N1()->DeadP() || (*w)->N2()->DeadP()) ||
 	 ((*w)->C1() == (*w)->C2() ) ){
       // kill wall
 #ifdef QDEBUG
       qDebug() << "Killing wall." << endl;
 #endif
       (*w)->Kill();
 #ifdef QDEBUG
       if ((*w)) {
 	qDebug() << "Wall " << (*w)->Index() << " says: c1 = "
 		 << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl;
+      }
 #endif
       (*w)=0;
     } else {
 #ifdef QDEBUG
       qDebug() << "Not killing wall." << endl;
       qDebug() << "Wall " << (*w)->Index() << " says: c1 = "
 	       << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl;
 #endif
+    }
+  }
   walls.remove(0);
   // Unregister me from my nodes, and delete the node if it no longer belongs to any cells
   list<Node *> superfluous_nodes;
   for (list<Node *>::iterator n=nodes.begin();
        n!=nodes.end();
        n++) {
   for (list<Node *>::iterator n=nodes.begin(); n!=nodes.end(); n++) {
     Node &no(*(*n));
     // locate myself in the node's owner list
     list<Neighbor>::iterator cellpos;
     bool cell_found=false;
     for (list<Neighbor>::iterator nb=no.owners.begin();
 	 nb!=no.owners.end();
 	 nb++) {
     for (list<Neighbor>::iterator nb=no.owners.begin(); nb!=no.owners.end(); nb++) {
       if (nb->cell == this) {
 	cellpos = nb;
 	cell_found = true;
 	break;
+      }
+    }
     if (!cell_found) {
       // I think this cannot happen; but if I am wrong, unpredictable things would happen. So throw an exception.
       throw ("Cell not found in CellBase::Apoptose()\n\rPlease correct the code to handle this situation.");
+    }
     Neighbor tmp = *cellpos;
     no.owners.erase(cellpos);
     // if node has no owners left, or only has a connection to special cell -1 (outside world), mark it as dead.
     if (no.owners.size()==0 || (no.owners.size()==1 && no.owners.front().cell->BoundaryPolP()) ) {
       no.MarkDead();
     } else {
       // register node with outside world
       if (find_if( no.owners.begin(), no.owners.end(),
 		   bind2nd ( mem_fun_ref(&Neighbor::CellEquals), m->boundary_polygon->Index() ) ) == no.owners.end() ) {
 	tmp.cell = m->boundary_polygon;
 	no.owners.push_back(tmp);
+      }
+    }
+  }
   // mark cell as dead
   MarkDead();
+}
 void Cell::ConstructConnections(void)
+{
   // Tie up the nodes of this cell, assuming they are correctly ordered
   //cerr << "Constructing connections of cell " << index << endl;
   for (list<Node *>::iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (list<Node *>::iterator i=nodes.begin(); i!=nodes.end(); i++) {
     //cerr << "Connecting node " << *i << endl;
     //cerr << "Node " << *i << endl << " = " << *(*i) << endl;
     // 1. Tidy up existing connections (which are part of this cell)
     if ((*i)->owners.size()>0) {
       list<Neighbor>::iterator neighb_with_this_cell=
 	// remove myself from the neighbor list of the node
 	find_if((*i)->owners.begin(),
 		(*i)->owners.end(),
 		bind2nd(mem_fun_ref( &Neighbor::CellEquals ),this->Index() )  );
       if (neighb_with_this_cell!=(*i)->owners.end())
 	(*i)->owners.erase(neighb_with_this_cell);
+    }
     Node *previous;
     if (i!=nodes.begin()) {
       list<Node *>::iterator previous_iterator=i;
       previous_iterator--;
       previous=*previous_iterator;
     } else {
       previous=nodes.back();
+    }
     Node *next;
     list<Node *>::iterator next_iterator=i;
     next_iterator++;
     if (next_iterator==nodes.end()) {
       next=nodes.front();
     } else {
       next=*next_iterator;
+    }
     (*i)->owners.push_back( Neighbor( this, previous, next ) );
+  }
+}
 bool Cell::DivideOverGivenLine(const Vector v1, const Vector v2, bool fix_cellwall, NodeSet *node_set )
+{
   if (dead) return false;
   // check each edge for intersection with the line
   ItList new_node_locations;
 #ifdef QDEBUG
   qDebug() << "Cell " << Index() << " is doing DivideOverGivenLine" << endl;
 #endif
   for (list<Node *>::iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (list<Node *>::iterator i=nodes.begin(); i!=nodes.end(); i++) {
     Vector v3 = *(*i);
     list<Node *>::iterator nb=i;
     nb++;
     if (nb == nodes.end()) {
       nb = nodes.begin();
+    }
     Vector v4 = *(*nb);
     double denominator =
       (v4.y - v3.y)*(v2.x - v1.x) - (v4.x - v3.x)*(v2.y - v1.y);
     double ua =
       ((v4.x - v3.x)*(v1.y - v3.y) - (v4.y - v3.y)*(v1.x -v3.x))/denominator;
     double ub =
       ((v2.x - v1.x)*(v1.y-v3.y) - (v2.y- v1.y)*(v1.x - v3.x))/denominator;
     //cerr << "Edge " << *i << " to " << *nb << ": ua = " << ua << ", ub = " << ub << ":  ";
     // this construction with "TINY" should simulate open/closed interval <0,1]
     if ( ( TINY < ua && ua < 1.+TINY ) && ( TINY < ub && ub < 1.+TINY ) ) {
       // yes, intersection detected. Push the location to the list of iterators
       new_node_locations.push_back(nb);
+    }
+}
+  }
 #ifdef QDEBUG
   if (new_node_locations.size()<2) {
     qDebug() << "Line does not intersect with two edges of Cell " << Index() << endl;
     qDebug() << "new_node_locations.size() = " << new_node_locations.size() << endl;
     return false;
+  }
 ItList::iterator i = new_node_locations.begin();
 list< Node *>::iterator j;
 qDebug() << "-------------------------------" << endl;
 qDebug() << "Location of new nodes: " << (**i)->Index() << " and ";
   ItList::iterator i = new_node_locations.begin();
   list< Node *>::iterator j;
   qDebug() << "-------------------------------" << endl;
   qDebug() << "Location of new nodes: " << (**i)->Index() << " and ";
 ++i;
 j = *i;
 if (j==nodes.begin()) j=nodes.end(); j--;
   ++i;
   j = *i;
   if (j==nodes.begin()) j=nodes.end(); j--;
 qDebug() << (*j)->Index() << endl;
 qDebug() << "-------------------------------" << endl;
   qDebug() << (*j)->Index() << endl;
   qDebug() << "-------------------------------" << endl;
 if ( **new_node_locations.begin() == *j ) {
   qDebug() << "Rejecting proposed division (cutting off zero area)." << endl;
   return false;
+ }
   if ( **new_node_locations.begin() == *j ) {
     qDebug() << "Rejecting proposed division (cutting off zero area)." << endl;
     return false;
+  }
 #endif
 DivideWalls(new_node_locations, v1, v2, fix_cellwall, node_set);
+  DivideWalls(new_node_locations, v1, v2, fix_cellwall, node_set);
 return true;
+  return true;
+}
 // Core division procedure
 void Cell::DivideWalls(ItList new_node_locations, const Vector from, const Vector to, bool fix_cellwall, NodeSet *node_set)
+{
   if (dead) return;
   bool boundary_touched_flag=false;
   // Step 0: keep some data about the parent before dividing
   ParentInfo parent_info;
   parent_info.polarization = ReduceCellAndWalls<Vector>( PINdir );
   parent_info.polarization.Normalise();
   parent_info.PINmembrane = SumTransporters(1);
   parent_info.PINendosome = Chemical(1);
   //cerr << "Parent polarization before division: " << parent_info.polarization << endl;
   // Step 1: create a daughter cell
   Cell *daughter=m->AddCell(new Cell());
   // Step 2: Copy the basics of parent cell to daughter
   for (int i=0;i<NChem();i++) {
     daughter->chem[i]=chem[i];
+  }
   daughter->cell_type = cell_type;
   for (int i=0;i<NChem();i++) {
     daughter->new_chem[i]=new_chem[i];
+  }
   daughter->boundary=boundary;
   daughter->m=m;
   daughter->target_area=target_area/2.;
   target_area/=2;
   daughter->cellvec=cellvec;
   // Division currently only works for convex cells: i.e. if the division line
   // intersects the cells at two points only.
   if (new_node_locations.size()!=2) {
     // Note: if you would add the possibility of dividing non-convex
     // cells, remember to update the code below. There are some
     // fixed-size arrays over there!
     cerr << "Warning in Cell::Division: division of non-convex cells not fully implemented" << endl;
     // Reject the daughter cell and decrement the amount of cells
     // again. We can do this here because it is the last cell added.
     // Never, ever try to fully erase a cell elsewhere, because we
     // make heavy use of cell indices in this project; if you erase a
     // Cell somewhere in the middle of Mesh::Cells the indices will
     // get totally messed up...! (e.g. the indices used in Nodes::cells)
 #ifdef QDEBUG
     qDebug() << "new_node_locations.size() = " << new_node_locations.size() <<endl;
     qDebug() << "daughter->index = " << daughter->index << endl;
     qDebug() << "cells.size() = " << m->cells.size() << endl;
 #endif
     m->cells.pop_back();
     Cell::NCells()--;
     m->shuffled_cells.pop_back();
     return;
+  }
   // We can be sure we only need two positions here because divisions
   // of non-convex cells are rejected above.
   Vector new_node[2];
   Node *new_node_ind[2];
   int new_node_flag[2];
   Edge div_edges[2];
   int nnc=0;
   Wall *div_wall[4];
   double orig_length[2];
   for (int i=0;i<4;i++) { div_wall[i]=0; orig_length[i/2] = 0.; }
   // construct new Nodes at the intersection points
   // unless they coincide with existing points
   for ( ItList::const_iterator i=new_node_locations.begin();
 	i!=new_node_locations.end();
 	i++) {
   for ( ItList::const_iterator i=new_node_locations.begin(); i!=new_node_locations.end(); i++) {
     // intersection between division axis
     // and line from this node to its predecessor
     // method used: http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
     Vector v1 = from;
     Vector v2 = to;
     Vector v3 = *(**i);
     // get previous node
     list<Node *>::iterator nb=*i;
     if (nb == nodes.begin()) {
       nb = nodes.end();
+    }
     nb--;
     Vector v4=*( *nb );
     double denominator =
       (v4.y - v3.y)*(v2.x - v1.x) - (v4.x - v3.x)*(v2.y - v1.y);
     double ua =
       ((v4.x - v3.x)*(v1.y - v3.y) - (v4.y - v3.y)*(v1.x -v3.x))/denominator;
     double intersec_x = v1.x + ua*(v2.x-v1.x);
     double intersec_y = v1.y + ua*(v2.y-v1.y);
     // construct a new node at intersec
     // we construct a vector temporarily,
     // until we are sure we are going to keep the node...
     // Node count is damaged if we construct superfluous nodes
     Vector *n=new Vector(intersec_x,intersec_y,0);
     div_edges[nnc].first=*nb;
     div_edges[nnc].second=**i;
     // Insert this new Node if it is far enough (5% of element length)
     // from one of the two existing nodes, else use existing node
     //
     // old, fixed value was: par.collapse_node_threshold = 0.05
     double collapse_node_threshold = 0.05;
 #ifdef FLEMING
     collapse_node_threshold = par.collapse_node_threshold;
 #endif
     double elem_length = ( (*(**i)) - (*(*nb)) ).Norm();
     if ( ( *(**i) - *n ).Norm() < collapse_node_threshold  * elem_length ) {
       new_node_flag[nnc]=1;
       new_node[nnc] = *(**i);
       new_node_ind[nnc] = **i;
       //cerr << **i << endl ;
     } else
       if ( (*(*nb) - *n).Norm() < collapse_node_threshold * elem_length ) {
 	new_node_flag[nnc]=2;
 	new_node[nnc] = *(*nb);
 	new_node_ind[nnc] = *nb;
       } else {
 	new_node_flag[nnc]=0;
 	new_node[nnc] = *n;
+      }
     nnc++;
     delete n;
+  }
   for (int i=0;i<2;i++) {
     Cell *neighbor_cell=0; // we need this to split up the "Wall" objects.
     // for both divided edges:
     //      insert its new node into all cells that own the divided edge
     // but only if it really is a new node:
     if (new_node_flag[i]!=0) {
       if (fix_cellwall) {
 	(new_node_ind[i])->fixed = true;
 	// all this we'll do later for the node set :-)
 	/* (new_node_ind[i])->boundary = true;
 	   (new_node_ind[i])->sam = true;
 	   boundary = SAM;
 	   daughter->boundary = SAM;
 	   boundary_touched_flag = true;
 	*/
+      }
     } else {
       // (Construct a list of all owners:)
       // really construct the new node (if this is a new node)
       new_node_ind[i] =
 	m->AddNode(new Node (new_node[i]) );
       // if a new node is inserted into a fixed edge (i.e. in the petiole)
       // make the new node fixed as well
@@ @@ -601,195 +587,193 @@ void Cell::DivideWalls(ItList new_node_l @@
 #ifdef QDEBUG
       list<Neighbor> unique_owners;
       copy(owners.begin(), owners.end(), back_inserter(unique_owners));
       unique_owners.unique( mem_fun_ref( &Neighbor::Eq ) );
       qDebug() << "The dividing edge nodes: " << div_edges[i].first->Index()
 	       << " and " << div_edges[i].second->Index() << " are owned by cells: ";
       // spit out each owners' cell index
       foreach(Neighbor neighbor, unique_owners){
 	qDebug() << neighbor.cell->Index() << "  ";
+      }
       qDebug() << endl;
 #endif
       // Search through the sorted list of edge node owners looking for duplicate pairs. Each pair represents an actual edge owner.
       list<Neighbor> edge_owners;
       list<Neighbor>::iterator it;
       for (it=owners.begin(); it!=owners.end(); it++) {
 	it = adjacent_find(it, owners.end(), neighbor_cell_eq);
 	if (it == owners.end()) break; // bail if reach the end of the list
 #ifdef QDEBUG
 	qDebug() << "Considering: " << it->cell->Index() << " as a possible edge owner." << endl;
 #endif
 	if (it->cell->Index() != this->Index()) {
 #ifdef QDEBUG
 	  qDebug() << "Adding: " << it->cell->Index() << " to the list of edge owners." << endl;
 #endif
 	  edge_owners.push_back(*it);
+	}
+      }
       if (edge_owners.size() > 1){
 	// Remove the boundary polygon - if its there
 	list<Neighbor>::iterator it;
 	if ((it = find_if (edge_owners.begin(), edge_owners.end(), bind2nd(mem_fun_ref(&Neighbor::CellEquals), -1)))
 	    != edge_owners.end()) {
 #ifdef QDEBUG
 	  qDebug() << "deleating: " << it->cell->Index() << " from the list of edge owners." << endl;
 #endif
 	  edge_owners.erase(it);
+	}
+      }
 #ifdef QDEBUG
       qDebug() << "The edge owners list has: " << edge_owners.size() << " elements" << endl;
 #endif
       // Since the list should always contain exactly one element, pass it on as an iterator
       list<Neighbor>::iterator c = (edge_owners.size() != 0) ? edge_owners.begin() : edge_owners.end();
       // (can we have more than one neighboring cell here??)
       if (c!=owners.end()) {
 	neighbor_cell = c->cell;
 	if (!c->cell->BoundaryPolP()) {
 	  // find correct position in the cells node list
 	  // to insert the new node
 	  list<Node *>::iterator ins_pos = find
 	    (neighbor_cell->nodes.begin(),
 	     neighbor_cell->nodes.end(),
 	     div_edges[i].first);
 	  neighbor_cell->nodes.insert(ins_pos, new_node_ind[i]);
 	  neighbor_cell->ConstructConnections();
 	  // give walls to daughter later
+	}
       } else {
 	neighbor_cell = 0;
+      }
+    }
     // Split the Wall with the neighboring cell
     // if the neighbor cell has not yet been identified above, do it now
     if (neighbor_cell == 0) {
       list<Neighbor> owners;
       // push all cells owning the two nodes of the divides edges
       // onto a list
       copy((div_edges[i].first)->owners.begin(),
 	   (div_edges[i].first)->owners.end(),
 	   back_inserter(owners));
       copy((div_edges[i].second)->owners.begin(),
 	   (div_edges[i].second)->owners.end(),
 	   back_inserter(owners));
       // find first non-self duplicate in the owners:
       // cells owning the same two nodes
       // share an edge with me
       owners.sort( mem_fun_ref( &Neighbor::Cmp ) );
       list<Neighbor>::iterator c;
       for (c=owners.begin();
 	   c!=owners.end();
 	   c++) {
       for (c=owners.begin(); c!=owners.end(); c++) {
 	c=adjacent_find(c,owners.end(),neighbor_cell_eq);
 	if (c->cell->Index() != this->Index() || c==owners.end()) break;
+      }
       if (c!=owners.end())
 	neighbor_cell = c->cell;
       else
 	neighbor_cell = 0;
+    }
     if (neighbor_cell /* && !neighbor_cell->BoundaryPolP() */) {
       //cerr << "Cell "  << index << " says: neighboring cell is " << neighbor_cell->index << endl;
       /*************** 1. Find the correct wall element  ********************/
       list<Wall *>::iterator w, start_search;
       w = start_search = walls.begin();
       do {
 	// Find wall between this cell and neighbor cell
 	w = find_if( start_search, walls.end(), bind2nd (mem_fun( &Wall::is_wall_of_cell_p ), neighbor_cell ) );
 	start_search = w; start_search++; // continue searching at next element
       } while ( w!=walls.end() && !(*w)->IntersectsWithDivisionPlaneP( from, to ) ); // go on until we find the right one.
       if (w == walls.end()) {
 #ifdef QDEBUG
 	qDebug() << "Whoops, wall element not found...!" << endl;
 	qDebug() << "Cell ID: " << neighbor_cell->Index() << endl;
 	qDebug() << "My cell ID: " << Index() << endl;
 #endif
       } else {
 	// 2. Split it up, if we should (sometimes, the new node coincides with an existing node so
 	// we should not split up the Wall)
 	if (new_node_ind[i]!=(*w)->n1 && new_node_ind[i]!=(*w)->n2) {
 	  Wall *new_wall;
 	  // keep the length of the original wall; we need it to equally divide the transporter concentrations
 	  // over the two daughter walls
 	  (*w)->SetLength(); // make sure we've got the current length
 	  orig_length[i] = (*w)->Length();
 	  //cerr << "Original length is " << orig_length[i] << endl;
 	  if ((*w)->c1 == this ) {
 	    //  cerr << "Cell " << (*w)->c1->Index() << " splits up wall " << *(*w) << ", into: " << endl;
 	    new_wall = new Wall( (*w)->n1, new_node_ind[i], this, neighbor_cell);
 	    (*w)->n1 = new_node_ind[i];
 	    //  cerr << "wall " << *(*w) << ", and new wall " << *new_wall << endl;
 	  } else {
 	    new_wall = new Wall( (*w)->n1, new_node_ind[i], neighbor_cell, this);
 	    (*w)->n1 = new_node_ind[i];
+	  }
 	  // 3. Give wall elements to appropriate cells
 	  if (new_wall->n1 != new_wall->n2) {
 	    if (par.copy_wall)
 	      new_wall->CopyWallContents(**w);
 	    else {
 	      // If wall contents are not copied, decide randomly which wall will be the "parent"
 	      // otherwise we will get biases (to the left), for example in the meristem growth model
 	      if (RANDOM()<0.5) {
 		new_wall->SwapWallContents(*w);
+	      }
+	    }
 	    AddWall(new_wall);
 	    // cerr << "Building new wall: this=" << Index() << ", neighbor_cell = " << neighbor_cell->Index() << endl;
 	    neighbor_cell->AddWall( new_wall);
 	    //cerr << "Existing wall: c1 = " << (*w)->c1->Index() << ", neighbor_cell = " << (*w)->c2->Index() << endl;
 	    // Remember the addresses of the new walls
 	    div_wall[2*i+0] = *w;
 	    div_wall[2*i+1] = new_wall;
 	    // we will correct the transporter concentrations later in this member function, after division
 	    // First the new nodes should be inserted into the cells, before we can calculate wall lengths
 	    // Remember that cell walls can be bent, so have a bigger length than the Euclidean distance n1->n2
 	  } else {
 	    delete new_wall;
+	  }
+	}
+      }
+    }
   }  // closing loop over the two divided edges (for (int i=0;i<2;i++) )
   // move half of the nodes to the daughter
+  {
@@ @@ -884,893 +868,820 @@ void Cell::DivideWalls(ItList new_node_l @@
 #ifdef MULTIPLE_NODES
   // intermediate, extra nodes
   // Calculate distance between the two new nodes
   double dist=( new_node[1] - new_node[0] ).Norm();
   //bool fixed_wall = (new_node_ind[0])->fixed && (new_node_ind[1])->fixed;
   bool fixed_wall = false;
   // Estimate number of extra nodes in wall
   // factor 4 is to keep tension on the walls;
   // this is a hidden parameter and should be made explicit
   // later on.
   int n=(int)((dist/Node::target_length)/4+0.5);
   Vector nodevec = ( new_node[1]- new_node[0]).Normalised();
   double element_length = dist/(double)(n+1);
   // note that wall nodes need to run in inverse order in parent
   list<Node *>::iterator ins_pos = daughter->nodes.end();
   for (int i=1;i<=n;i++) {
     Node *node=
       m->AddNode( new Node( new_node[0] + i*element_length*nodevec ) );
     node->fixed=fixed_wall;
     if (!fix_cellwall)
       node->boundary = false;
     else {
       node->fixed = true;
+    }
     ins_pos=daughter->nodes.insert(ins_pos, node );
     new_nodes_parent.push_back( node );
     // optionally add the new node to the nodeset (passed by pointer)
     // (in this way we can move the NodeSet as a whole; useful for a fixed cutting line)
     if (node_set) {
       node_set->AddNode( node );
+    }
+  }
 #endif
   daughter->nodes.push_back( new_node_ind[0] );
   new_nodes_parent.push_back( new_node_ind[1] );
   // optionally add the new node to the nodeset (passed by pointer)
   // (in this way we can move the NodeSet as a whole; useful for a fixed cutting line)
   if (node_set) {
     node_set->AddNode( new_node_ind[1] );
+  }
   // move the new nodes to the parent
   nodes.clear();
   copy( new_nodes_parent.begin(),
 	new_nodes_parent.end(),
 	back_inserter(nodes) );
   // Repair cell lists of Nodes, and node connectivities
   ConstructConnections();
   daughter->ConstructConnections();
   if (boundary_touched_flag) {
     m->boundary_polygon->ConstructConnections();
+  }
   // collecting neighbors of divided cell
   list<CellBase *> broken_neighbors;
   // this cell's old neighbors
   copy(neighbors.begin(), neighbors.end(), back_inserter(broken_neighbors) );
   // this cell
   broken_neighbors.push_back(this);
   // its daughter
   broken_neighbors.push_back(daughter);
   // Recalculate area of parent and daughter
   area = CalcArea();
   daughter->area = daughter->CalcArea();
   SetIntegrals();
   daughter->SetIntegrals();
   // Add a "Cell Wall" for diffusion algorithms
   Wall *wall = new Wall( new_node_ind[0], new_node_ind[1], this, daughter );
   AddWall( wall );
   daughter->AddWall( wall );
   //cerr << "Correct walls of cell " << Index() << " and daughter " << daughter->Index() << endl;
   // Move Walls to daughter cell
   list <Wall *> copy_walls = walls;
   for (list<Wall *>::iterator w = copy_walls.begin();
        w!=copy_walls.end();
        w++) {
   for (list<Wall *>::iterator w = copy_walls.begin(); w!=copy_walls.end(); w++) {
     //cerr << "Doing wall, before:  " << **w << endl;
     //  checks the nodes of the wall and gives it away if appropriate
     (*w)->CorrectWall ( );
     //cerr << "and after: " << **w << endl;
+  }
   // Correct tranporterconcentrations of divided walls
   for (int i=0;i<4;i++) {
     if (div_wall[i]) {
       div_wall[i]->SetLength();
       div_wall[i]->CorrectTransporters(orig_length[i/2]);
+    }
+  }
   //cerr << "Cell " << index << " has been dividing, and gave birth to Cell " << daughter->index << endl;
   // now reconstruct neighbor list for all "broken" neighbors
   for (list<CellBase *>::iterator i=broken_neighbors.begin();
        i!=broken_neighbors.end();i++) {
   for (list<CellBase *>::iterator i=broken_neighbors.begin(); i!=broken_neighbors.end(); i++) {
     ((Cell *)(*i))->ConstructNeighborList();
+  }
   ConstructNeighborList();
   daughter->ConstructNeighborList();
   m->plugin->OnDivide(&parent_info, daughter, this);
   daughter->div_counter=(++div_counter);
+}
 // Move the whole cell
 void Cell::Move(const Vector T) {
   for (list<Node *>::const_iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (list<Node *>::const_iterator i=nodes.begin(); i!=nodes.end(); i++) {
     *(*i)+=T;
+  }
+}
 double Cell::Displace(double dx, double dy, double dh)
+{
   // Displace whole cell, add resulting energy to dh,
   // and accept displacement if energetically favorable
   //
   // Method is called if a "fixed" node is displaced
   // Warning: length constraint not yet  CORRECTLY implemented for this function
   // Attempt to move this cell in a random direction
   //  Vector movement(par.mc_cell_stepsize*(RANDOM()-0.5),par.mc_cell_stepsize*(RANDOM()-0.5),0);
   dh=0;
   Vector movement(dx,dy,0);
   vector< pair<Node *, Node *> > length_edges;
   vector<double> cellareas;
   cellareas.reserve(neighbors.size());
   // for the length constraint, collect all edges to this cell's nodes,
   // which are not part of the cell
   // the length of these edges will change
   double old_length=0.;
   for (list<Node *>::const_iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
     for (list<Neighbor>::const_iterator n=(*i)->owners.begin();
 	 n!=(*i)->owners.end();
 	 n++) {
   for (list<Node *>::const_iterator i=nodes.begin(); i!=nodes.end(); i++) {
     for (list<Neighbor>::const_iterator n=(*i)->owners.begin(); n!=(*i)->owners.end(); n++) {
       if (n->getCell()!=this) {
 	length_edges.push_back( pair <Node *,Node *> (*i, n->nb1) );
 	length_edges.push_back( pair <Node *,Node *> (*i, n->nb2) );
 	old_length +=
 	  DSQR(Node::target_length-(*(*i)-*(n->nb1)).Norm())+
 	  DSQR(Node::target_length-(*(*i)-*(n->nb2)).Norm());
+      }
+    }
+  }
   // calculate area energy difference of neighboring cells
   // (this cells' shape remains unchanged)
   double old_area_energy=0., old_length_energy=0.;
   for (list<CellBase *>::const_iterator i=neighbors.begin();
        i!=neighbors.end();
        i++) {
   for (list<CellBase *>::const_iterator i=neighbors.begin(); i!=neighbors.end(); i++) {
     old_area_energy += DSQR((*i)->Area()-(*i)->TargetArea());
     old_length_energy += DSQR((*i)->Length()-(*i)->TargetLength());
+  }
   Move(movement);
   double new_area_energy=0., new_length_energy=0.;
   for (list<CellBase *>::const_iterator i=neighbors.begin();
        i!=neighbors.end();
        i++) {
   for (list<CellBase *>::const_iterator i=neighbors.begin(); i!=neighbors.end(); i++) {
     cellareas.push_back((*i)->CalcArea());
     new_area_energy += DSQR(cellareas.back()-(*i)->TargetArea());
     new_length_energy += DSQR((*i)->CalcLength()-(*i)->TargetLength());
+  }
   double new_length=0;
   for ( vector< pair< Node *, Node * > >::const_iterator e = length_edges.begin();
 	e != length_edges.end();
 	e++) {
     new_length +=  DSQR(Node::target_length-
 			(*(e->first)-*(e->second)).Norm());
   for ( vector< pair< Node *, Node * > >::const_iterator e = length_edges.begin(); e != length_edges.end(); e++) {
     new_length +=  DSQR(Node::target_length - (*(e->first) - *(e->second)).Norm());
+  }
   dh += (new_area_energy - old_area_energy) + (new_length_energy - old_length_energy) * lambda_celllength +
     par.lambda_length * (new_length - old_length);
   if (dh<0 || RANDOM()<exp(-dh/par.T)) {
     // update areas of cells
     //cerr << "neighbors: ";
     list<CellBase *>::const_iterator nb_it = neighbors.begin();
     for (vector<double>::const_iterator ar_it = cellareas.begin();
 	 ar_it!=cellareas.end();
 	 ( ar_it++, nb_it++) ) {
     for (vector<double>::const_iterator ar_it = cellareas.begin(); ar_it!=cellareas.end(); ( ar_it++, nb_it++) ) {
       ((Cell *)(*nb_it))->area = *ar_it;
       (*nb_it)->SetIntegrals();
+    }
     //cerr << endl;
   } else {
     Move ( -1*movement);
+  }
   return dh;
+}
 void Cell::Displace (void)
+{
   Displace(par.mc_cell_stepsize*(RANDOM()-0.5),par.mc_cell_stepsize*(RANDOM()-0.5),0);
+}
 // Get energy level of whole cell (excluding length constraint?)
 double Cell::Energy(void) const
+{
   double energy = 0.;
   double length_contribution = 0.;
   for (list<Node *>::const_iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
     for (list<Neighbor>::const_iterator n=(*i)->owners.begin();
 	 n!=(*i)->owners.end();
 	 n++) {
   for (list<Node *>::const_iterator i=nodes.begin(); i!=nodes.end(); i++) {
     for (list<Neighbor>::const_iterator n=(*i)->owners.begin(); n!=(*i)->owners.end(); n++) {
       if (n->getCell()==this) {
 	length_contribution +=
 	  DSQR(Node::target_length-(*(*i)-*(n->nb1)).Norm())+
+	  DSQR(Node::target_length-(*(*i)-*(n->nb1)).Norm()) +
 	  DSQR(Node::target_length-(*(*i)-*(n->nb2)).Norm());
+      }
+    }
+  }
   // wall elasticity constraint
   energy += par.lambda_length * length_contribution;
   // area constraint
   energy += DSQR(CalcArea() - target_area);
   // cell length constraint
   energy += lambda_celllength * DSQR(Length() - target_length);
   return energy;
+}
 bool Cell::SelfIntersect(void)
+{
   // The (obvious) O(N*N) algorithm
   // Compare each edge against each other edge
   // An O(N log(N)) algorithm by Shamos & Hoey (1976) supposedly exists;
   // it was mentioned on comp.graphics.algorithms
   // But I haven't been able to lay my hand on the paper.
   // Let's try whether we need it....
   // method used: http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
   for (list<Node *>::const_iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (list<Node *>::const_iterator i=nodes.begin(); i!=nodes.end(); i++) {
     list<Node *>::const_iterator j=i;
     ++j;
     for (;
 	 j!=nodes.end();
 	 j++)
     for (; j!=nodes.end(); j++)
+      {
 	Vector v1 = *(*i);
 	list<Node *>::const_iterator nb=i;
 	nb++;
 	if (nb == nodes.end()) {
 	  nb = nodes.begin();
+	}
 	Vector v2 = *(*nb);
 	Vector v3 = *(*j);
 	nb=j;
 	nb++;
 	if (nb == nodes.end()) {
 	  nb = nodes.begin();
+	}
 	Vector v4=*( *nb );
 	double denominator =
 	  (v4.y - v3.y)*(v2.x - v1.x) - (v4.x - v3.x)*(v2.y - v1.y);
 	double ua =
 	  ((v4.x - v3.x)*(v1.y - v3.y) - (v4.y - v3.y)*(v1.x -v3.x))/denominator;
 	double ub =
 	  ((v2.x - v1.x)*(v1.y-v3.y) - (v2.y- v1.y)*(v1.x - v3.x))/denominator;
 	if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) {
 	  //cerr << "ua = " << ua << ", ub = " << ub << endl;
 	  return true;
+	}
+      }
+  }
   return false;
+}
 bool Cell::MoveSelfIntersectsP(Node *moving_node_ind, Vector new_pos)
+{
   // Check whether the polygon will self-intersect if moving_node_ind
   // were displaced to new_pos
   // Compare the two new edges against each other edge
   // O(2*N)
   // method used for segment intersection:
   // http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
   Vector neighbor_of_moving_node[2];
   //cerr << "list<Node *>::const_iterator moving_node_ind_pos = find (nodes.begin(),nodes.end(),moving_node_ind);\n";
   list<Node *>::const_iterator moving_node_ind_pos = find (nodes.begin(),nodes.end(),moving_node_ind);
   list<Node *>::const_iterator nb = moving_node_ind_pos;
   //cerr << "Done\n";
   nb++;
   if (nb == nodes.end()) {
     nb = nodes.begin();
+  }
   neighbor_of_moving_node[0]=*(*nb);
   nb=moving_node_ind_pos;
   if (nb == nodes.begin()) {
     nb = nodes.end();
+  }
   nb--;
   neighbor_of_moving_node[1]=*( *nb );
   for (list<Node *>::const_iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (list<Node *>::const_iterator i=nodes.begin(); i!=nodes.end(); i++) {
     for (int j=0;j<2;j++) { // loop over the two neighbors of moving node
       list<Node *>::const_iterator nb=i;
       nb++;
       if (nb == nodes.end()) {
 	nb = nodes.begin();
+      }
       if (*i == moving_node_ind || *nb == moving_node_ind) {
 	// do not compare to self
 	continue;
+      }
       Vector v3 = *(*i);
       Vector v4 = *(*nb);
       double denominator =
 	(v4.y - v3.y)*(neighbor_of_moving_node[j].x - new_pos.x) - (v4.x - v3.x)*(neighbor_of_moving_node[j].y - new_pos.y);
       double ua =
 	((v4.x - v3.x)*(new_pos.y - v3.y) - (v4.y - v3.y)*(new_pos.x -v3.x))/denominator;
       double ub =
 	((neighbor_of_moving_node[j].x - new_pos.x)*(new_pos.y-v3.y) - (neighbor_of_moving_node[j].y- new_pos.y)*(new_pos.x - v3.x))/denominator;
       if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) {
 	//cerr << "ua = " << ua << ", ub = " << ub << endl;
 	return true;
+      }
+    }
+  }
   return false;
+}
 /*! \brief Test if this cell intersects with the given line.
  */
 bool Cell::IntersectsWithLineP(const Vector v1, const Vector v2)
+{
   // Compare the line against each edge
   // method used: http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
   for (list<Node *>::const_iterator i=nodes.begin(); i!=nodes.end(); i++)
+    {
       Vector v3 = *(*i);
       list<Node *>::const_iterator nb=i;
       nb++;
       if (nb == nodes.end()) {
 	nb = nodes.begin();
+      }
       Vector v4 = *(*nb);
       double denominator =
 	(v4.y - v3.y)*(v2.x - v1.x) - (v4.x - v3.x)*(v2.y - v1.y);
       double ua =
 	((v4.x - v3.x)*(v1.y - v3.y) - (v4.y - v3.y)*(v1.x -v3.x))/denominator;
       double ub =
 	((v2.x - v1.x)*(v1.y-v3.y) - (v2.y- v1.y)*(v1.x - v3.x))/denominator;
       if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) {
 	return true;
+      }
+    }
   return false;
+}
 /*! \brief Constructs Walls, but only one per cell boundary.
   Standard method constructs a Wall for each cell wall element,
   making transport algorithms computationally more intensive than needed.
   We can remove this? Well, let's leave it in the code in case we need it for something else. E.g. for importing leaf architectures in different formats than our own... :-)
 */
 void Cell::ConstructWalls(void)
+{
   return;
   if (dead) return;
   walls.clear();
   neighbors.clear();
   // Get "corner points; i.e. nodes where more than 2 cells are connected
   list<Node *> corner_points;
   for (list<Node *>::const_iterator i=nodes.begin();
        i!=nodes.end();i++) {
   for (list<Node *>::const_iterator i=nodes.begin(); i!=nodes.end();i++) {
     // look for nodes belonging to >2 cells
     if ((*i)->owners.size()>2) {
       // push onto list
       corner_points.push_back(*i);
+    }
+  }
   // Construct Walls between corner points
   // previous one in list
   list<Node *>::const_iterator nb = (--corner_points.end());
   // loop over list,
   for (list<Node *>::const_iterator i=corner_points.begin();
        i!=corner_points.end(); ( i++, nb++) ) {
   for (list<Node *>::const_iterator i=corner_points.begin(); i!=corner_points.end(); ( i++, nb++) ) {
     if (nb==corner_points.end()) nb=corner_points.begin();
     // add owning cells to a list
     list<Cell *> owning_cells;
     Node &n(*(*i));
     for (list<Neighbor>::const_iterator j=n.owners.begin();
 	 j!=n.owners.end();
 	 j++) {
     for (list<Neighbor>::const_iterator j=n.owners.begin(); j!=n.owners.end(); j++) {
       owning_cells.push_back(j->cell);
+    }
     Node &n2(*(*nb));
     for (list<Neighbor>::const_iterator j=n2.owners.begin();
 	 j!=n2.owners.end();
 	 j++) {
     for (list<Neighbor>::const_iterator j=n2.owners.begin(); j!=n2.owners.end(); j++) {
       owning_cells.push_back(j->cell);
+    }
     // sort cell owners
     owning_cells.sort( mem_fun( &Cell::Cmp ));
     // find duplicates
     vector<Cell *> duplicates;
     list<Cell *>::const_iterator prevj = (--owning_cells.end());
     for (list<Cell *>::const_iterator j=owning_cells.begin();
 	 j!=owning_cells.end();
 	 ( j++, prevj++) ) {
       if (prevj==owning_cells.end()) prevj=owning_cells.begin();
       if (*j==*prevj) duplicates.push_back(*j);
     for (list<Cell *>::const_iterator j=owning_cells.begin(); j!=owning_cells.end(); ( j++, prevj++) ) {
       if (prevj==owning_cells.end())
 	prevj=owning_cells.begin();
       if (*j==*prevj)
 	duplicates.push_back(*j);
+    }
     if (duplicates.size()==3) { // ignore cell boundary (this occurs only after the first division, I think)
       vector<Cell *>::iterator dup_it = find_if(duplicates.begin(),duplicates.end(),mem_fun(&Cell::BoundaryPolP) );
       if (dup_it!=duplicates.end())
 	duplicates.erase(dup_it);
       else {
 	return;
+      }
+    }
     // One Wall for each neighbor, so we should be able to correctly construct neighbor lists here.
     if (duplicates[0]==this) {
       AddWall(  new Wall(*nb,*i,duplicates[0],duplicates[1]) );
       if (!duplicates[1]->BoundaryPolP()) {
 	neighbors.push_back(duplicates[1]);
+      }
     } else {
       AddWall ( new Wall(*nb,*i,duplicates[1],duplicates[0]) );
       if (!duplicates[0]->BoundaryPolP()) {
 	neighbors.push_back(duplicates[0]);
+      }
+    }
+  }
+}
 void BoundaryPolygon::Draw(QGraphicsScene *c, QString tooltip)
+{
   // Draw the BoundaryPolygon on a QCanvas object
   CellItem* p = new CellItem(this, c);
   QPolygonF pa(nodes.size());
   int cc=0;
   for (list<Node *>::const_iterator n=nodes.begin();
        n!=nodes.end();
        n++) {
   for (list<Node *>::const_iterator n=nodes.begin(); n!=nodes.end(); n++) {
     Node *i=*n;
     pa[cc++] = QPoint((int)((Offset().x+i->x)*Factor()),
 		      (int)((Offset().y+i->y)*Factor()) );
     pa[cc++] = QPoint((int)((Offset().x+i->x)*Factor()), (int)((Offset().y+i->y)*Factor()) );
+  }
   p->setPolygon(pa);
   p->setPen(par.outlinewidth>=0?QPen( QColor(par.cell_outline_color),par.outlinewidth):QPen(Qt::NoPen));
   p->setBrush( Qt::NoBrush );
   p->setZValue(1);
   if (!tooltip.isEmpty())
     p->setToolTip(tooltip);
   p->show();
+}
 void Cell::Flux(double *flux, double *D)
+{
   // loop over cell edges
   for (int c=0;c<NChem();c++) flux[c]=0.;
   for (int c=0;c<NChem();c++)
     flux[c]=0.;
   for (list<Wall *>::iterator i=walls.begin();
        i!=walls.end();
        i++) {
   for (list<Wall *>::iterator i=walls.begin(); i!=walls.end(); i++) {
     // leaf cannot take up chemicals from environment ("no flux boundary")
     if ((*i)->c2->BoundaryPolP()) continue;
     // flux depends on edge length and concentration difference
     for (int c=0;c<NChem();c++) {
       double phi = (*i)->length * ( D[c] ) * ( ((Cell *)(*i)->c2)->chem[c] - chem[c] );
 #ifdef QDEBUG
       if ((*i)->c1!=this) {
 	qDebug() << "Warning, bad cells boundary: " << (*i)->c1->Index() << ", " << index << endl;
+      }
 #endif
       flux[c] += phi;
+    }
+  }
+}
 // graphics stuff, not compiled for batch versions
 #ifdef QTGRAPHICS
 #include "canvas.h"
 void Cell::Draw(QGraphicsScene *c, QString tooltip)
+{
   // Draw the cell on a QCanvas object
   if (DeadP()) {
 #ifdef QDEBUG
     qDebug() << "Cell " << index << " not drawn, because dead." << endl;
 #endif
     return;
+  }
   CellItem* p = new CellItem(this, c);
   QPolygonF pa(nodes.size());
   int cc=0;
   for (list<Node *>::const_iterator n=nodes.begin();
        n!=nodes.end();
        n++) {
   for (list<Node *>::const_iterator n=nodes.begin(); n!=nodes.end(); n++) {
     Node *i=*n;
     pa[cc++] = QPoint((int)((offset[0]+i->x)*factor),
 		      (int)((offset[1]+i->y)*factor) );
+  }
   QColor cell_color;
   m->plugin->SetCellColor(this,&cell_color);
   p->setPolygon(pa);
   p->setPen(par.outlinewidth>=0?QPen( QColor(par.cell_outline_color),par.outlinewidth):QPen(Qt::NoPen));
   p->setBrush( cell_color );
   p->setZValue(1);
   if (!tooltip.isEmpty())
     p->setToolTip(tooltip);
   p->show();
+}
 void Cell::DrawCenter(QGraphicsScene *c) const {
   // Maginfication derived similarly to that in nodeitem.cpp
   // Why not use Cell::Magnification()?
   const double mag = par.node_mag;
   // construct an ellipse
   QGraphicsEllipseItem *disk = new QGraphicsEllipseItem ( -1*mag, -1*mag, 2*mag, 2*mag, 0, c );
   disk->setBrush( QColor("forest green") );
   disk->setZValue(5);
   disk->show();
   Vector centroid=Centroid();
   disk -> setPos((offset[0]+centroid.x)*factor,(offset[1]+centroid.y)*factor);
+}
 void Cell::DrawNodes(QGraphicsScene *c) const {
   for (list<Node *>::const_iterator n=nodes.begin();
        n!=nodes.end();
        n++) {
   for (list<Node *>::const_iterator n=nodes.begin(); n!=nodes.end(); n++) {
     Node *i=*n;
     NodeItem *item = new NodeItem ( &(*i), c );
     item->setColor();
     item->setZValue(5);
     item->show();
     item ->setPos(((offset[0]+i->x)*factor),
 		  ((offset[1]+i->y)*factor) );
     item ->setPos(((offset[0]+i->x)*factor), ((offset[1]+i->y)*factor) );
+  }
+}
 void Cell::DrawIndex(QGraphicsScene *c) const {
   DrawText( c, QString("%1").arg(index));
+}
 // Draw any text in the cell's center
 void Cell::DrawText(QGraphicsScene *c, const QString &text) const {
   Vector centroid = Centroid();
   QGraphicsSimpleTextItem *ctext = new QGraphicsSimpleTextItem ( text, 0, c );
   ctext->setPen( QPen(QColor(par.textcolor)) );
   ctext->setZValue(20);
   ctext->setFont( QFont( "Helvetica", par.cellnumsize, QFont::Bold) );
   ctext->show();
   ctext ->setPos(((offset[0]+centroid.x)*factor),
 		 ((offset[1]+centroid.y)*factor) );
+}
 void Cell::DrawAxis(QGraphicsScene *c) const {
   Vector long_axis;
   double width;
   Length(&long_axis, &width);
   //cerr << "Length is "  << length << endl;
   long_axis.Normalise();
   Vector short_axis=long_axis.Perp2D();
   Vector centroid = Centroid();
   Vector from = centroid - 0.5 * width * short_axis;
   Vector to = centroid + 0.5 * width *short_axis;
   QGraphicsLineItem *line = new QGraphicsLineItem(0, c);
   line->setPen( QPen(QColor(par.arrowcolor),2) );
   line->setZValue(2);
   line->setLine( ( (offset[0]+from.x)*factor ),
 		 ( (offset[1]+from.y)*factor ),
 		 ( (offset[0]+to.x)*factor ),
 		 ( (offset[1]+to.y)*factor ) );
   line->setZValue(10);
   line->show();
+}
 void Cell::DrawStrain(QGraphicsScene *c) const {
   c = NULL; // assignment merely to obviate compilation warning
   MyWarning::warning("Sorry, Cell::DrawStrain temporarily not implemented.");
+}
 void Cell::DrawFluxes(QGraphicsScene *c, double arrowsize)
+{
   // get the mean flux through this cell
   Vector vec_flux = ReduceCellAndWalls<Vector>( PINdir );
   vec_flux.Normalise();
   vec_flux *= arrowsize;
   QGraphicsArrowItem *arrow = new QGraphicsArrowItem(0,c);
   Vector centroid = Centroid();
   Vector from = centroid - vec_flux/2.;
   Vector to = centroid + vec_flux/2.;
   arrow->setPen( QPen(QColor(par.arrowcolor),par.outlinewidth));
   arrow->setZValue(2);
   arrow->setLine( ( (offset[0]+from.x)*factor ),
 		  ( (offset[1]+from.y)*factor ),
 		  ( (offset[0]+to.x)*factor ),
 		  ( (offset[1]+to.y)*factor ) );
   arrow->setZValue(10);
   arrow->show();
+}
 void Cell::DrawWalls(QGraphicsScene *c) const {
   for_each(walls.begin(), walls.end(), bind2nd ( mem_fun ( &Wall::Draw ) , c ) );
   // to see the cells connected the each wall (for debugging), uncomment the following
   //for_each(walls.begin(), walls.end(), bind2nd ( mem_fun ( &Wall::ShowStructure ), c ) );
+}
 void Cell::DrawValence(QGraphicsScene *c) const {
   DrawText(c, QString("%1").arg(walls.size()) );
+}
 #endif // QTGRAPHICS !
 /*! \brief Recalculate the lengths of the cell's Walls.
   Call this function after the Monte Carlo updates, and before doing the reaction-diffusion iterations.
 */
 void Cell::SetWallLengths(void)
+{
   for (list<Wall *>::iterator de=walls.begin();
        de!=walls.end();
        de++) {
   for (list<Wall *>::iterator de=walls.begin(); de!=walls.end(); de++) {
     // Step 1: find the path of nodes leading along the Wall.
     // A Wall often represents a curved cell wall: we want the total
     // length _along_ the wall here...
     // Locate first and second nodes of the edge in list of nodes
     list<Node *>::const_iterator first_node_edge = find(nodes.begin(), nodes.end(), (*de)->n1);
     list<Node *>::const_iterator second_node_edge_plus_1 = ++find(nodes.begin(), nodes.end(), (*de)->n2);
     double sum_length = 0.;
     // Now, walk to the second node of the edge in the list of nodes
     for (list<Node *>::const_iterator n=++first_node_edge;
 	 n!=second_node_edge_plus_1;
 	 ++n ) {
       if (n==nodes.end()) n=nodes.begin(); /* wrap around */
     for (list<Node *>::const_iterator n=++first_node_edge; n!=second_node_edge_plus_1; ++n ) {
       if (n==nodes.end())
 	n=nodes.begin(); /* wrap around */
       list<Node *>::const_iterator prev_n = n;
       if (prev_n==nodes.begin()) prev_n=nodes.end();
       if (prev_n==nodes.begin())
 	prev_n=nodes.end();
       --prev_n;
       // Note that Node derives from a Vector, so we can do vector calculus as defined in vector.h
       sum_length += (*(*prev_n) - *(*n)).Norm();
       //cerr << "Node " << *prev_n << " to " << *n << ", cumulative length = " << sum_length << endl;
+    }
     // We got the total length of the Wall now, store it:
     (*de)->length = sum_length;
     //cerr << endl;
     // goto next de
+  }
+}
 //! Add Wall w to the list of Walls
 void Cell::AddWall( Wall *w )
+{
   // if necessary, we could try later inserting it at the correct position
 #ifdef QDEBUG
   if (w->c1 == w->c2 ){
     qDebug() << "Wall between identical cells: " << w->c1->Index()<< endl;
+  }
 #endif
   // Add Wall to Cell's list
   walls.push_back( w );
   // Add wall to Mesh's list if it isn't there yet
   if (find (
 	    m->walls.begin(), m->walls.end(),
 	    w )
       == m->walls.end() ) {
     m->walls.push_back(w);
   if (find ( m->walls.begin(), m->walls.end(), w ) == m->walls.end() ) { m->walls.push_back(w);
+  }
+}
 //! Remove Wall w from the list of Walls
 list<Wall *>::iterator Cell::RemoveWall( Wall *w )
+{
   // remove wall from Mesh's list
   m->walls.erase(
 		 find(
 		      m->walls.begin(), m->walls.end(),
 		      w )
 		 );
   m->walls.erase( find( m->walls.begin(), m->walls.end(), w ) );
   // remove wall from Cell's list
   return walls.erase (find( walls.begin(), walls.end(), w ));
+}
 void Cell::EmitValues(double t)
+{
   //  cerr << "Attempting to emit " << t << ", " << chem[0] << ", " << chem[1] << endl;
   emit ChemMonValue(t, chem);
+}
 /* finis */

src/cellbase.cpp

➞

Show inline comments

@@ @@ -141,510 +141,489 @@ CellBase::CellBase(const CellBase &src) @@
   chem=new double[NChem()];
   for (int i=0;i<NChem();i++) {
     chem[i]=src.chem[i];
+  }
   new_chem=new double[NChem()];
   for (int i=0;i<NChem();i++) {
     new_chem[i]=src.new_chem[i];
+  }
   boundary=src.boundary;
   area=src.area;
   target_length=src.target_length;
   lambda_celllength=src.lambda_celllength;
   intgrl_xx=src.intgrl_xx; intgrl_xy=src.intgrl_xy; intgrl_yy=src.intgrl_yy;
   intgrl_x=src.intgrl_x; intgrl_y=src.intgrl_y;
   target_area=src.target_area;
   index=src.index;
   nodes=src.nodes;
   neighbors=src.neighbors;
   walls=src.walls;
   source = src.source;
   fixed = src.fixed;
   source_conc = src.source_conc;
   source_chem = src.source_chem;
   cellvec = src.cellvec;
   at_boundary=src.at_boundary;
   pin_fixed = src.pin_fixed;
   stiffness = src.stiffness;
   marked = src.marked;
   dead = src.dead;
   cell_type = src.cell_type;
   div_counter = src.div_counter;
   flag_for_divide = src.flag_for_divide;
   division_axis = src.division_axis;
+}
 CellBase CellBase::operator=(const CellBase &src)
+{
   Vector::operator=(src);
   for (int i=0;i<NChem();i++) {
     chem[i]=src.chem[i];
+  }
   for (int i=0;i<NChem();i++) {
     new_chem[i]=src.chem[i];
+  }
   boundary=src.boundary;
   area=src.area;
   intgrl_xx=src.intgrl_xx; intgrl_xy=src.intgrl_xy; intgrl_yy=src.intgrl_yy;
   intgrl_x=src.intgrl_x; intgrl_y=src.intgrl_y;
   target_area=src.target_area;
   target_length=src.target_length;
   lambda_celllength=src.lambda_celllength;
   index=src.index;
   nodes=src.nodes;
   neighbors=src.neighbors;
   walls=src.walls;
   source = src.source;
   fixed = src.fixed;
   source_conc = src.source_conc;
   source_chem = src.source_chem;
   cellvec = src.cellvec;
   at_boundary=src.at_boundary;
   pin_fixed = src.pin_fixed;
   stiffness = src.stiffness;
   marked = src.marked;
   dead = src.dead;
   cell_type = src.cell_type;
   div_counter = src.div_counter;
   flag_for_divide = src.flag_for_divide;
   division_axis = src.division_axis;
   return *this;
+}
 void CellBase::SetChemical(int c, double conc)
+{
   if (c>=NChem()) {
     stringstream error;
     error << "SetChemical: value c = " << c << " is out of range\n";
     throw error.str().c_str();
+  }
   chem[c]=conc;
+}
 void CellBase::SetTransporters(int ch, double conc)
+{
   if (ch>=NChem()) {
     stringstream error;
     error << "SetChemical: value ch = " << ch << " is out of range\n";
     throw error.str().c_str();
+  }
   for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
   for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
     (*w)->setTransporter(this, ch, conc);
+  }
+}
 ostream &CellBase::print(ostream &os) const
+{
   os << "[ index = " << index << " {" << x << ", " << y << ", " << z << "}: {";
   for (int i=0;i<NChem()-1;i++) {
     os << chem[i] << ", ";
+  }
   os << chem[NChem()-1] << " } ]";
   os << endl << "Nodelist = { " << endl;
   for (list<Node *>::const_iterator i =  nodes.begin(); i!=nodes.end(); i++) {
     os << (*i)->Index() << "( " << *i << ") ";
+  }
   os << " } ";
   for (list<Wall *>::const_iterator i =  walls.begin(); i!=walls.end(); i++) {
     (*i)->print(os);
     os << ", ";
+  }
   os << endl;
   os << " [ area = " << area << " ]";
   os << " [ walls = ";
   for (list<Wall *>::const_iterator i= walls.begin();
        i!=walls.end();
        i++) {
   for (list<Wall *>::const_iterator i= walls.begin(); i!=walls.end(); i++) {
     os << (*i)->n1->Index() << " -> " << (*i)->n2->Index() << ", " <<  (*i)->c1->Index() << " | " << (*i)->c2->Index() << ", ";
+  }
   os << " ] ";
   os << "div_counter = " << div_counter << endl;
   os << "cell_type = " << cell_type << endl;
   os << endl;
   return os;
+}
 ostream &operator<<(ostream &os, const CellBase &c)
+{
   c.print(os);
   return os;
+}
 double CellBase::CalcArea(void) const
+{
   double loc_area=0.;
   for (list<Node *>::const_iterator i=nodes.begin();
        i!=(nodes.end());
        i++) {
   for (list<Node *>::const_iterator i=nodes.begin(); i!=(nodes.end()); i++) {
     list<Node *>::const_iterator i_plus_1=i; i_plus_1++;
     if (i_plus_1==nodes.end())
       i_plus_1=nodes.begin();
     loc_area+= (*i)->x * (*i_plus_1)->y;
     loc_area-= (*i_plus_1)->x * (*i)->y;
+  }
   // http://technology.niagarac.on.ca/courses/ctec1335/docs/arrays2.pdf
   return fabs(loc_area)/2.0;
+}
 Vector CellBase::Centroid(void) const
+{
   double area=0.;
   double integral_x_dxdy=0.,integral_y_dxdy=0.;
   for (list<Node *>::const_iterator i=nodes.begin();
        i!=(nodes.end());
        i++) {
   for (list<Node *>::const_iterator i=nodes.begin(); i!=(nodes.end()); i++) {
     list<Node *>::const_iterator i_plus_1=i; i_plus_1++;
     if (i_plus_1==nodes.end())
       i_plus_1=nodes.begin();
     area+= (*i)->x * (*i_plus_1)->y;
     area-= (*i_plus_1)->x * (*i)->y;
     integral_x_dxdy+=
       ((*i_plus_1)->x+(*i)->x)*
       ((*i)->x*(*i_plus_1)->y-
        (*i_plus_1)->x*(*i)->y);
     integral_y_dxdy+=
       ((*i_plus_1)->y+(*i)->y)*
       ((*i)->x*(*i_plus_1)->y-
        (*i_plus_1)->x*(*i)->y);
+  }
   area = fabs(area)/2.0;
   integral_x_dxdy/=6.;
   integral_y_dxdy/=6.;
   Vector centroid(integral_x_dxdy,integral_y_dxdy,0);
   centroid/=area;
   return centroid;
+}
 void CellBase::SetIntegrals(void) const
+{
   // Set the initial values for the integrals over x^2,
   // xy, yy, x, and y
   // these values will be updated after each move of the CellBase wall
   intgrl_xx=0.; intgrl_xy=0.; intgrl_yy=0.;
   intgrl_x=0.; intgrl_y=0.;
   area=0.;
   list<Node *>::const_iterator nb;
   list<Node *>::const_iterator i=nodes.begin();
   for (;
        i!=(nodes.end());
        i++) {
   for (; i!=(nodes.end()); i++) {
     nb = i; nb++; if (nb==nodes.end()) nb=nodes.begin();
     area+=(*i)->x*(*nb)->y;
     area-=(*nb)->x*(*i)->y;
     intgrl_xx+=
       ((*i)->x*(*i)->x+
        (*nb)->x*(*i)->x+
        (*nb)->x*(*nb)->x ) *
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y);
     intgrl_xy+=
       ((*nb)->x*(*i)->y-
        (*i)->x*(*nb)->y)*
       ((*i)->x*(2*(*i)->y+(*nb)->y)+
        (*nb)->x*((*i)->y+2*(*nb)->y));
     intgrl_yy+=
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y)*
       ((*i)->y*(*i)->y+
        (*nb)->y*(*i)->y+
        (*nb)->y*(*nb)->y );
     intgrl_x+=
       ((*nb)->x+(*i)->x)*
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y);
     intgrl_y+=
       ((*nb)->y+(*i)->y)*
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y);
+  }
   area = fabs(area)/2.0;
+}
 double CellBase::Length(Vector *long_axis, double *width)  const
+{
   // Calculate length and axes of CellBase
   // Calculate inertia tensor
   // see file inertiatensor.nb for explanation of this method
   if (!lambda_celllength) {
     // Without length constraint we do not keep track of the cells'
     // moments of inertia. So we must calculate them here.
     SetIntegrals();
+  }
   double intrx=intgrl_x/6.;
   double intry=intgrl_y/6.;
   double ixx=(intgrl_xx/12.)-(intrx*intrx)/area;
   double ixy=(intgrl_xy/24.)+(intrx*intry)/area;
   double iyy=(intgrl_yy/12.)-(intry*intry)/area;
   double rhs1=(ixx+iyy)/2., rhs2=sqrt( (ixx-iyy)*(ixx-iyy)+4*ixy*ixy )/2.;
   double lambda_b=rhs1+rhs2;
   // see: http://scienceworld.wolfram.com/physics/MomentofInertiaEllipse.html
   //    cerr << "n = " << n << "\n";
   if (long_axis) {
     *long_axis = Vector(-ixy, lambda_b - ixx, 0);
     //   cerr << "ixx = " << ixx << ", ixy = " << ixy << ", iyy = " << iyy << ", area = " << area << endl;
+  }
   if (width) {
     *width = 4*sqrt((rhs1-rhs2)/area);
+  }
   return 4*sqrt(lambda_b/area);
+}
 double CellBase::CalcLength(Vector *long_axis, double *width)  const
+{
   // Calculate length and axes of CellBase, without touching cells raw moments
   // Calculate inertia tensor
   // see file inertiatensor.nb for explanation of this method
   double my_intgrl_xx=0., my_intgrl_xy=0., my_intgrl_yy=0.;
   double my_intgrl_x=0., my_intgrl_y=0., my_area=0.;
   my_area=0.;
   list<Node *>::const_iterator nb;
   list<Node *>::const_iterator i=nodes.begin();
   for (;
        i!=(nodes.end());
        i++) {
   for (; i!=(nodes.end()); i++) {
     nb = i; nb++; if (nb==nodes.end()) nb=nodes.begin();
     my_area+=(*i)->x*(*nb)->y;
     my_area-=(*nb)->x*(*i)->y;
     my_intgrl_xx+=
       ((*i)->x*(*i)->x+
        (*nb)->x*(*i)->x+
        (*nb)->x*(*nb)->x ) *
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y);
     my_intgrl_xy+=
       ((*nb)->x*(*i)->y-
        (*i)->x*(*nb)->y)*
       ((*i)->x*(2*(*i)->y+(*nb)->y)+
        (*nb)->x*((*i)->y+2*(*nb)->y));
     my_intgrl_yy+=
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y)*
       ((*i)->y*(*i)->y+
        (*nb)->y*(*i)->y+
        (*nb)->y*(*nb)->y );
     my_intgrl_x+=
       ((*nb)->x+(*i)->x)*
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y);
     my_intgrl_y+=
       ((*nb)->y+(*i)->y)*
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y);
+  }
   //my_area/=2.0;
   my_area = fabs(my_area)/2.0;
   double intrx=my_intgrl_x/6.;
   double intry=my_intgrl_y/6.;
   double ixx=(my_intgrl_xx/12.)-(intrx*intrx)/my_area;
   double ixy=(my_intgrl_xy/24.)+(intrx*intry)/my_area;
   double iyy=(my_intgrl_yy/12.)-(intry*intry)/my_area;
   double rhs1=(ixx+iyy)/2., rhs2=sqrt( (ixx-iyy)*(ixx-iyy)+4*ixy*ixy )/2.;
   double lambda_b=rhs1+rhs2;
   // see: http://scienceworld.wolfram.com/physics/MomentofInertiaEllipse.html
   //    cerr << "n = " << n << "\n";
   if (long_axis) {
     *long_axis = Vector(-ixy, lambda_b - ixx, 0);
     //   cerr << "ixx = " << ixx << ", ixy = " << ixy << ", iyy = " << iyy << ", my_area = " << my_area << endl;
+  }
   if (width) {
     *width = 4*sqrt((rhs1-rhs2)/my_area);
+  }
   return 4*sqrt(lambda_b/my_area);
+}
 void CellBase::ConstructNeighborList(void)
+{
   neighbors.clear();
   for (//list<Wall *>::const_reverse_iterator wit=walls.rbegin();
        list<Wall *>::const_iterator wit=walls.begin();
        // somehow the reverse_iterator returns by walls needs to be casted to const to let this work.
        // it seems to me it is a bug in the STL implementation...
        wit!=walls.end();
        wit++) {
     if ((*wit)->C1() != this) {
       neighbors.push_back((*wit)->C1());
     } else {
       neighbors.push_back((*wit)->C2());
+    }
+  }
   // remove all boundary_polygons from the list
   list <CellBase *>::iterator e=neighbors.begin();
   at_boundary=false;
   do {
     // Code crashes here after cutting off part of the leaf. I can't find the problem.
     // Leaving the "Illegal" walls in the simulation helps. (c1=-1 && c2=-1)
     // Work-around: define leaf primordium. Save to XML. Restart. Read XML file.
     // Sorry about this; I hope to solve this annoying issue later. RM :-).
     // All cells in neighbors seem to be okay (I might be messing some part of the memory elsewhere
     // during the cutting operation?).
     e = find_if(neighbors.begin(),neighbors.end(),mem_fun(&CellBase::BoundaryPolP));
     if (e!=neighbors.end()) {
       e=neighbors.erase(e);
       at_boundary=true;
     } else {
       break;
+    }
   } while(1);
+}
 // Save the cell to a stream so we can reconstruct its state later
 void CellBase::Dump(ostream &os) const
+{
   os << index << " " << nodes.size() << endl;
   Vector::Dump(os);
   os << endl;
   for (list<Node *>::const_iterator i=nodes.begin();i!=nodes.end();i++) {
     os << *i << " ";
+  }
   os << endl;
   os << index << " " << neighbors.size() << endl;
   for (list<CellBase *>::const_iterator i=neighbors.begin();i!=neighbors.end();i++) {
     os << *i << " ";
+  }
   os << endl << walls.size() << endl << endl;
   os << NChem() << " ";
   for (int i=0;i<NChem();i++) {
     os << chem[i] << " ";
+  }
   os << endl;
   os << NChem() << " ";
   for (int i=0;i<NChem();i++) {
     os << new_chem[i] << " ";
+  }
   os << endl;
   os << boundary << " " << area << " " << target_area << " " << target_length
      << " " << fixed << " " << intgrl_xx << " " << intgrl_xy << " " << intgrl_yy
      << " " << intgrl_x << " " << intgrl_y << " " << source << " ";
   cellvec.Dump(os);
   os << " " << source_conc << " " << source_chem;
   os << endl;
+}
 void CellBase::UnfixNodes(void)
+{
   for (list<Node *>::const_iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (list<Node *>::const_iterator i=nodes.begin(); i!=nodes.end(); i++) {
     (*i)->Unfix();
+  }
+}
 void CellBase::FixNodes(void)
+{
   for (list<Node *>::const_iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (list<Node *>::const_iterator i=nodes.begin(); i!=nodes.end(); i++) {
     (*i)->Fix();
+  }
+}
 // returns true if cell is at border
 bool CellBase::AtBoundaryP(void) const
+{
   return at_boundary;
+}
 QString CellBase::printednodelist(void)
+{
   QString info_string = "Nodelist = { ";
   for (list<Node *>::const_iterator i =  nodes.begin(); i!=nodes.end(); i++) {
     info_string += QString("%1 ").arg((*i)->Index());
+  }
   info_string += " } ";
   return info_string;
+}
 /* finis*/

src/mesh.cpp

➞

Show inline comments

 /*
+ *
  *  This file is part of the Virtual Leaf.
+ *
  *  The Virtual Leaf is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
+ *
  *  The Virtual Leaf is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU General Public License for more details.
+ *
  *  You should have received a copy of the GNU General Public License
  *  along with the Virtual Leaf.  If not, see <http://www.gnu.org/licenses/>.
+ *
  *  Copyright 2010 Roeland Merks.
+ *
  */
 #include <string>
 #include <algorithm>
 #include <vector>
 #include <sstream>
 #include <cstdlib>
 //#include <cerrno>
 #include <cstring>
 #include <numeric>
 #include <functional>
 #include <fstream>
 #include <QPair>
 #include "mesh.h"
 #include "tiny.h"
 #include "parameter.h"
 #include "random.h"
 #include "pi.h"
 #include "parse.h"
 #include "matrix.h"
 #include "sqr.h"
 #include "nodeset.h"
 #include "nodeitem.h"
 #include "simplugin.h"
 #include <QDebug>
 #include <set>
 #include <iostream>
 #include <iterator>
 static const std::string _module_id("$Id$");
 extern Parameter par;
 void Mesh::AddNodeToCellAtIndex(Cell *c, Node *n, Node *nb1, Node *nb2, list<Node *>::iterator ins_pos) {
   c->nodes.insert(ins_pos, n);
   n->owners.push_back( Neighbor(c, nb1, nb2 ) );
+}
 void Mesh::AddNodeToCell(Cell *c, Node *n, Node *nb1, Node *nb2) {
   c->nodes.push_back( n );
   n->owners.push_back( Neighbor(c, nb1, nb2 ) );
+}
 void Mesh::PerturbChem(int chemnum, double range) {
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
   for (vector<Cell *>::iterator i=cells.begin(); i!=cells.end(); i++) {
     (*i)->chem[chemnum] += range*(RANDOM()-0.5);
     if ((*i)->chem[chemnum]<0.) (*i)->chem[chemnum]=0.;
+  }
+}
 void Mesh::CellFiles(const Vector ll, const Vector ur) {
   Cell *cell = RectangularCell(ll,ur,0.001);
   for (int c=0;c<Cell::NChem();c++) {
     cell->SetChemical(c,par.initval[c]);
+  }
   cell->SetTargetArea(cell->CalcArea());
   Vector axis(1,0,0);
   // divide rectangle a number of times
   for (int i=0;i<6;i++) {
     IncreaseCellCapacityIfNecessary();
     vector <Cell *> current_cells = cells;
     for (vector<Cell *>::iterator j=current_cells.begin();
 	 j!=current_cells.end();j++) {
     for (vector<Cell *>::iterator j=current_cells.begin(); j!=current_cells.end();j++) {
       (*j)->DivideOverAxis(axis);
+    }
     axis=axis.Perp2D();
+  }
   IncreaseCellCapacityIfNecessary();
   axis=axis.Perp2D();
   vector <Cell *> current_cells = cells;
   for (vector<Cell *>::iterator j=current_cells.begin();
        j!=current_cells.end();j++) {
   for (vector<Cell *>::iterator j=current_cells.begin(); j!=current_cells.end();j++) {
     (*j)->DivideOverAxis(axis);
+  }
   double sum_l=0; int n_l=0;
   for (list<Node *>::const_iterator i=cell->nodes.begin();
        i!=cell->nodes.end();
        i++) {
   for (list<Node *>::const_iterator i=cell->nodes.begin(); i!=cell->nodes.end(); i++) {
     list<Node *>::const_iterator nb=i; nb++;
     if (nb==cell->nodes.end())
       nb=cell->nodes.begin();
     double l = (**nb-**i).Norm();
     sum_l += l;
     n_l++;
+  }
   Node::target_length = sum_l/(double)n_l;
   // a bit more tension
   Node::target_length/=4.;
   SetBaseArea();
+}
 Cell *Mesh::RectangularCell(const Vector ll, const Vector ur, double rotation) {
   Cell *cell=AddCell(new Cell());
   cell->m=this;
   Matrix rotmat;
   rotmat.Rot2D(rotation); // rotation over 0,0
   Node *n1=AddNode(new Node(rotmat * ll));
   Node *n2=AddNode(new Node(rotmat * Vector(ll.x, ur.y,0)));
   Node *n3=AddNode(new Node(rotmat * ur));
   Node *n4=AddNode(new Node(rotmat * Vector(ur.x, ll.y,0)));
   n1->boundary=true;
   n2->boundary=true;
   n3->boundary=true;
   n4->boundary=true;
   //n1.fixed=n2.fixed=n3.fixed=n4.fixed=true;
   AddNodeToCell(cell, n4,
 		n1,
 		n3);
   AddNodeToCell(cell, n3,
 		n4,
 		n2);
   AddNodeToCell(cell, n2,
 		n3,
 		n1);
   AddNodeToCell(cell, n1,
 		n2,
 		n4);
   AddNodeToCell(boundary_polygon, n4,
 		n1,
 		n3);
   AddNodeToCell(boundary_polygon, n3,
 		n4,
 		n2);
   AddNodeToCell(boundary_polygon, n2,
 		n3,
 		n1);
   AddNodeToCell(boundary_polygon, n1,
 		n2,
 		n4);
   cell->setCellVec(Vector(0,1,0));
   boundary_polygon->m = this;
   boundary_polygon->area = 0;
   cell->area = cell->CalcArea();
   // target length is the length of the elements
   Node::target_length = ur.y-ll.y;
   // a bit more tension
   Node::target_length/=2;
   cell->SetIntegrals();
   cell->ConstructNeighborList();
   return cell;
+}
 Cell &Mesh::EllipticCell(double xc, double yc, double ra, double rb,  int nnodes, double rotation) {
   int first_node=Node::nnodes;
   //  nodes.reserve(nodes.size()+nnodes);
   //cells.push_back(Cell(xc,yc));
   Cell *c=AddCell(new Cell(xc,yc));
   c->m=this;
   for (int i=0;i<nnodes;i++) {
     double angle=2*Pi*(i/(double)nnodes);
     double x=xc+ra*cos(angle)*cos(rotation) - rb*sin(angle)*sin(rotation);
     double y=yc+ra*cos(angle)*sin(rotation) + rb*sin(angle)*cos(rotation);
     Node *n=AddNode(new Node(x,y,0));
     n->boundary = true;
+  }
   for (int i=0;i<nnodes;i++) {
     AddNodeToCell(c,
 		  nodes[first_node + i],
 		  nodes[first_node+ (nnodes+i-1)%nnodes],
 		  nodes[first_node+ (i + 1)%nnodes]);
     AddNodeToCell(boundary_polygon,
 		  nodes[first_node + i],
 		  nodes[first_node+ (nnodes+i-1)%nnodes],
 		  nodes[first_node+ (i + 1)%nnodes]);
+  }
   boundary_polygon->m = this;
   boundary_polygon->area = 0;
   c->area = c->CalcArea();
   // target length is the length of the elements
   Node::target_length = (2 * ((ra +rb)/2.) * sin (Pi/nnodes));
   // a bit more tension
   Node::target_length/=2;
   //boundary_polygon = c;
   /*  list<int>::iterator nb;
       for (list<int>::iterator i=c->nodes.begin();
       i!=c->nodes.end();
       i++) {
       nb = i; nb++;
       if (nb==c->nodes.end()) {
       nb=c->nodes.begin();
+      }
       int next = *nb;
       nb = i;
       if (nb==c->nodes.begin()) {
       nb=c->nodes.end();
+      }
       nb--;
       int previous = *nb;
       getNode(*i).cells.push_back( Neighbor(boundary_polygon->index, next, previous) );
       }*/
   c->SetIntegrals();
   //c->ConstructNeighborList();
   //c->ConstructWalls();
   // initial cell has one wall with the outside world
   //c->walls.push_back( new Wall ( nodes.front(), nodes.front(), c, boundary_polygon ));
   c->at_boundary=true;
   return *c;
+}
 Cell &Mesh::LeafPrimordium(int nnodes, double pet_length) {
   // first leaf cell
   int first_node=Node::nnodes;
   Cell *circle=AddCell(new Cell(0,0));
   circle->m=this;
   const double ra=10, rb=10;
   const double xc=0,yc=0;
   const double rotation=0;
   for (int i=0;i<nnodes;i++) {
     double angle=2*Pi*(i/(double)nnodes);
     double x=xc+ra*cos(angle)*cos(rotation) - rb*sin(angle)*sin(rotation);
     double y=yc+ra*cos(angle)*sin(rotation) + rb*sin(angle)*cos(rotation);
     Node *n=AddNode(new Node(x,y,0));
     /* if (angle > 1.25*Pi && angle < 1.75*Pi ) {
        n.sam = true;
        }*/
     AddNodeToCell(circle,
 		  n,
 		  nodes[first_node+ (nnodes+i-1)%nnodes],
 		  nodes[first_node+ (i + 1)%nnodes]);
+  }
   boundary_polygon->m = this;
   boundary_polygon->area = 0;
   circle->area = circle->CalcArea();
   // target length is the length of the elements
   Node::target_length = (2 * ((ra +rb)/2.) * sin (Pi/nnodes));
   // a bit more tension
   Node::target_length/=2;
   circle->SetIntegrals();
   //return c;
   circle->SetTargetArea(2*circle->Area());
   // Petiole: starts at both sides of the circular cell
   // get position of the (n/4)'th and (3*(n/4))'th node.
   list<Node *>::reverse_iterator it_n1=circle->nodes.rbegin();
   for (int i=0;i<nnodes/2;i++)
+  for (int i=0; i<nnodes/2; i++)
     it_n1++;
   it_n1--;
   list<Node *>::reverse_iterator it_n2=--circle->nodes.rend();
   /* for (int i=0;i<n/2;i++)
      it_n2++;*/
   Cell *petiole=AddCell(new Cell());
   Node *n1 = *it_n1;
   Node *n2 = *it_n2;
   Node *n3=AddNode( new Node ( *n2 + Vector( 0, pet_length, 0) ) );
   Node *n4=AddNode( new Node ( *n1 + Vector( 0, pet_length, 0) ) );
   n3->boundary=true;
   n4->boundary=true;
   AddNodeToCell(petiole, *it_n1,
 		n4,
 		nodes[(*it_n2)->Index()
 		      + (( (*it_n1)->Index() - (*it_n2)->Index() )-1+nnodes)%nnodes]);
   list<Node *>::reverse_iterator i=it_n1; i++;
   for (;
        i!=it_n2;
        //(++i) == circle->nodes.rend() ? i : i=circle->nodes.rbegin() ) {
        ++i) {
     AddNodeToCell(petiole,
 		  *i,
   for (; i!=it_n2; ++i) {
     AddNodeToCell(petiole, *i,
 		  nodes[(*it_n2)->Index() + (((*i)->Index()-(*it_n2)->Index()) + 1)%nnodes],
 		  nodes[(*it_n2)->Index() + (((*i)->Index()-(*it_n2)->Index())-1+nnodes)%nnodes]);
+  }
   AddNodeToCell(petiole, *it_n2,
 		*it_n2 + 1,
 		n3);
   AddNodeToCell(petiole, *it_n2, *it_n2 + 1, n3);
   (*it_n2)->boundary=true;
   //petiole.nodes.push_back(n3.Index());
   //petiole.nodes.push_back(n4.Index());
   AddNodeToCell(petiole,
 		n3,
 		n2,
 		n4);
   AddNodeToCell(petiole,
 		n4,
 		n3,
 		n1);
   AddNodeToCell(petiole, n3, n2, n4);
   AddNodeToCell(petiole, n4, n3, n1);
 #ifdef QDEBUG
   qDebug() << circle << endl;
   qDebug() << petiole << endl;
 #endif
   AddNodeToCell(boundary_polygon, *it_n1,
 		n4,
 		*it_n2 + ((*it_n1-*it_n2)+1)%nnodes); // is this gonna work?
   AddNodeToCell(boundary_polygon, *it_n1, n4, *it_n2 + ((*it_n1-*it_n2)+1)%nnodes); // is this gonna work?
   (*it_n1)->boundary=true;
   for (int i=0;i<nnodes;i++) {
     if (nodes[(first_node + i)]->owners.size()==1) {
       AddNodeToCell(boundary_polygon,
 		    nodes[first_node +i],
 		    nodes[first_node+ (nnodes+i-1)%nnodes],
 		    nodes[first_node+ (i + 1)%nnodes]);
       nodes[first_node+i]->boundary=true;
+    }
+  }
   AddNodeToCell(boundary_polygon, *it_n2,
 		nodes[(nnodes+(*it_n2)->Index() - 1)%nnodes],
 		n3);
   AddNodeToCell(boundary_polygon,
 		n3,
 		n2,
 		n4);
   AddNodeToCell(boundary_polygon,
 		n4,
 		n3,
 		n1);
   AddNodeToCell(boundary_polygon, *it_n2, nodes[(nnodes+(*it_n2)->Index() - 1)%nnodes], n3);
   AddNodeToCell(boundary_polygon, n3, n2, n4);
   AddNodeToCell(boundary_polygon, n4, n3, n1);
   // make petiole solid
   for (list<Node *>::iterator i=petiole->nodes.begin();
        i!=petiole->nodes.end();
        i++) {
   for (list<Node *>::iterator i=petiole->nodes.begin(); i!=petiole->nodes.end(); i++) {
     (*i)->Fix();
+  }
   petiole->Fix();
   petiole->area=petiole->CalcArea();
   petiole->target_area=petiole->area;
   petiole->ConstructNeighborList();
   circle->ConstructNeighborList();
   boundary_polygon->ConstructConnections();
   boundary_polygon->ConstructNeighborList();
   circle->setCellVec(Vector(0,1,0));
   return *circle;
+}
 /*Cell &Mesh::Box() {
   }*/
 // return bounding box of mesh
 void Mesh::BoundingBox(Vector &LowerLeft, Vector &UpperRight) {
   LowerLeft = **nodes.begin();
   UpperRight = **nodes.begin();
   for (vector<Node *>::iterator c=nodes.begin();
        c!=nodes.end();
        c++) {
   for (vector<Node *>::iterator c=nodes.begin(); c!=nodes.end(); c++) {
     if ((*c)->x < LowerLeft.x)
       LowerLeft.x = (*c)->x;
     if ((*c)->y < LowerLeft.y)
       LowerLeft.y = (*c)->y;
     if ((*c)->z < LowerLeft.z)
       LowerLeft.z = (*c)->z;
     if ((*c)->x > UpperRight.x)
       UpperRight.x = (*c)->x;
     if ((*c)->y > UpperRight.y)
       UpperRight.y = (*c)->y;
     if ((*c)->z > UpperRight.z)
       UpperRight.z = (*c)->z;
+  }
+}
 double Mesh::Area(void) {
   double area=0;
   vector<Cell *>::iterator i=cells.begin();
   while (i != cells.end()) {
     area += (*(i++))->Area();
+  }
   return area;
+}
 void Mesh::SetBaseArea(void) {
   // Set base area to mean area.
   // This method is typically called during initiation, after
   // defining the first cell
   Cell::BaseArea() = Area()/cells.size();
+}
 // for optimization, we moved Displace to Mesh
 class DeltaIntgrl {
 public:
   double area;
   double ix, iy;
   double ixx,ixy,iyy;
   DeltaIntgrl(double sarea,double six,double siy,double sixx,double sixy,double siyy) {
     area=sarea;
     ix=six;
     iy=siy;
     ixx=sixx;
     ixy=sixy;
     iyy=siyy;
+  }
 };
 void Mesh::Clear(void) {
   // clear nodes
   for (vector<Node *>::iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (vector<Node *>::iterator i=nodes.begin(); i!=nodes.end(); i++) {
     delete *i;
+  }
   nodes.clear();
   Node::nnodes=0;
   node_insertion_queue.clear();
   // Clear NodeSets
   for (vector<NodeSet *>::iterator i=node_sets.begin();
        i!=node_sets.end();
        i++) {
   for (vector<NodeSet *>::iterator i=node_sets.begin(); i!=node_sets.end(); i++) {
     delete *i;
+  }
   node_sets.clear();
   time = 0;
   // clear cells
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
   for (vector<Cell *>::iterator i=cells.begin(); i!=cells.end(); i++) {
     delete *i;
+  }
   cells.clear();
   Cell::NCells() = 0;
   delete boundary_polygon;
   // Clear walls
   for (list<Wall *>::iterator i=walls.begin();
        i!=walls.end();
        i++) {
   for (list<Wall *>::iterator i=walls.begin(); i!=walls.end(); i++) {
     delete *i;
+  }
   walls.clear();
   WallBase::nwalls = 0;
   //tmp_walls->clear();
   shuffled_cells.clear();
   shuffled_nodes.clear();
   //Cell::ncells=0;
   /*    Cell::ClearNCells();
 	Node::nnodes=0;
 	cells.clear();
 	nodes.clear();
 	shuffled_cells.clear();
 	shuffled_nodes.clear();
 	node_insertion_queue.empty();
 	cerr << "Meshed cleared: cells: " << cells.size() << ", nodes: " << nodes.size() << endl;
   */
 #ifdef QDEBUG
   qDebug() << "cells.size() = " << cells.size() << endl;
   qDebug << "walls.size() = " << walls.size() << endl;
   qDebug << "nodes.size() = " << nodes.size() << endl;
 #endif
+}
 double Mesh::DisplaceNodes(void) {
   MyUrand r(shuffled_nodes.size());
   random_shuffle(shuffled_nodes.begin(),shuffled_nodes.end(),r);
   double sum_dh=0;
   list<DeltaIntgrl> delta_intgrl_list;
   for_each( node_sets.begin(), node_sets.end(), mem_fun( &NodeSet::ResetDone ) );
   for (vector<Node *>::const_iterator i=shuffled_nodes.begin();
        i!=shuffled_nodes.end();
        i++) {
   for (vector<Node *>::const_iterator i=shuffled_nodes.begin(); i!=shuffled_nodes.end(); i++) {
     //int n=shuffled_nodes[*i];
     Node &node(**i);
     // Do not allow displacement if fixed
     //if (node.fixed) continue;
     if (node.DeadP()) continue;
     // Attempt to move this cell in a random direction
     double rx=par.mc_stepsize*(RANDOM()-0.5); // was 100.
     double ry=par.mc_stepsize*(RANDOM()-0.5);
     // Uniform with a circle of radius par.mc_stepsize
     /* double r = RANDOM() * par.mc_stepsize;
        double th = RANDOM()*2*Pi;
        double rx = r * cos(th);
        double ry = r * sin(th);
     */
     Vector new_p(node.x+rx,node.y+ry,0);
     Vector old_p(node.x,node.y,0);
     /* if (node.boundary  && boundary_polygon->MoveSelfIntersectsP(n,  new_p )) {
     // reject if move of boundary results in self intersection
     continue;
     }*/
     if (node.node_set) {
       // move each node set only once
       if (!node.node_set->DoneP())
 	node.node_set->AttemptMove(rx,ry);
     } else {
       // for all cells to which this node belongs:
       //   calculate energy difference
       double area_dh=0.;
       double length_dh=0.;
       double bending_dh=0.;
       double cell_length_dh=0.;
       double alignment_dh=0.;
       double old_l1=0.,old_l2=0.,new_l1=0.,new_l2=0.;
       double sum_stiff=0.;
       double dh=0.;
       for (list<Neighbor>::const_iterator cit=node.owners.begin();
 	   cit!=node.owners.end();
 	   cit++) {
       for (list<Neighbor>::const_iterator cit=node.owners.begin(); cit!=node.owners.end(); cit++) {
 	//Cell &c=m->getCell(cit->cell);
 	//Cell &c=cit->cell->BoundaryPolP()?*boundary_polygon:*(cit->cell);
 	Cell &c=*((Cell *)(cit->cell));
 	//Cell &c=cells[cit->cell];
 	if (c.MoveSelfIntersectsP(&node,  new_p )) {
 	  // reject if move results in self intersection
 	  //
 	  // I know: using goto's is bad practice... except when jumping out
 	  // of deeply nested loops :-)
 	  //cerr << "Rejecting due to self-intersection\n";
 	  goto next_node;
+	}
 	// summing stiffnesses of cells. Move has to overcome this minimum required energy.
 	sum_stiff += c.stiffness;
 	// area - (area after displacement): see notes for derivation
 	//Vector i_min_1 = m->getNode(cit->nb1);
 	Vector i_min_1 = *(cit->nb1);
 	//Vector i_plus_1 = m->getNode(cit->nb2);
 	Vector i_plus_1 = *(cit->nb2);
 	// We must double the weights for the perimeter (otherwise they start bulging...)
 	double w1, w2;
 	if (node.boundary && cit->nb1->boundary)
 #ifdef FLEMING
 	  w1 = par.rel_perimeter_stiffness;
 #else
 	w1=2;
 #endif
 	else
 	  w1 = 1;
 	if (node.boundary && cit->nb2->boundary)
 #ifdef FLEMING
 	  w2 = par.rel_perimeter_stiffness;
 #else
 	w2 = 2;
 #endif
 	else
 	  w2 = 1;
 	//if (cit->cell>=0) {
 	if (!cit->cell->BoundaryPolP()) {
 	  double delta_A = 0.5 * ( ( new_p.x - old_p.x ) * (i_min_1.y - i_plus_1.y) +
 				   ( new_p.y - old_p.y ) * ( i_plus_1.x - i_min_1.x ) );
 	  area_dh +=  delta_A * (2 * c.target_area - 2 * c.area + delta_A);
 	  // cell length constraint
 	  // expensive and not always needed
 	  // so we check the value of lambda_celllength
 	  if (/* par.lambda_celllength */  cit->cell->lambda_celllength) {
 	    double delta_ix =
 	      (i_min_1.x + new_p.x)
 	      * (new_p.x * i_min_1.y- i_min_1.x * new_p.y) +
 	      (new_p.x + i_plus_1.x)
 	      * (i_plus_1.x * new_p.y- new_p.x * i_plus_1.y) -
 	      (i_min_1.x + old_p.x)
 	      * (old_p.x * i_min_1.y- i_min_1.x * old_p.y) -
 	      (old_p.x + i_plus_1.x)
 	      * (i_plus_1.x * old_p.y - old_p.x * i_plus_1.y);
 	    double delta_iy =
 	      (i_min_1.y + new_p.y)
 	      * (new_p.x * i_min_1.y- i_min_1.x * new_p.y) +
 	      (new_p.y + i_plus_1.y)
 	      * (i_plus_1.x * new_p.y- new_p.x * i_plus_1.y) -
 	      (i_min_1.y + old_p.y)
 	      * (old_p.x * i_min_1.y- i_min_1.x * old_p.y) -
 	      (old_p.y + i_plus_1.y)
 	      * (i_plus_1.x * old_p.y - old_p.x * i_plus_1.y);
 	    double delta_ixx =
 	      (new_p.x*new_p.x+
 	       i_min_1.x*new_p.x+
 	       i_min_1.x*i_min_1.x ) *
 	      (new_p.x*i_min_1.y - i_min_1.x*new_p.y) +
 	      (i_plus_1.x*i_plus_1.x+
 	       new_p.x*i_plus_1.x+
 	       new_p.x*new_p.x ) *
 	      (i_plus_1.x*new_p.y - new_p.x*i_plus_1.y) -
 	      (old_p.x*old_p.x+
 	       i_min_1.x*old_p.x+
 	       i_min_1.x*i_min_1.x ) *
@@ @@ -867,211 +785,205 @@ double Mesh::DisplaceNodes(void) { @@
 	     w1*(DSQR(new_l1)
 	     - DSQR(old_l1))
 	     + w2*(DSQR(new_l2)
 	     - DSQR(old_l2)); */
 	  /*	if (c.cellvec.ManhattanNorm()!=0) {
 	  // wall element length constraint
 	  double old_aniso1, old_aniso2;
 	  double new_aniso1, new_aniso2;
 	  // anisotropic expansion?
 	  old_aniso1 = 1 + par.lambda_aniso*(1 - fabs(InnerProduct((old_p-i_min_1), c.cellvec))/old_l1);
 	  old_aniso2 = 1 + par.lambda_aniso*(1 - fabs(InnerProduct((old_p-i_plus_1), c.cellvec))/old_l2);
 	  new_aniso1 = 1 + par.lambda_aniso*(1 - fabs(InnerProduct((new_p-i_min_1),  c.cellvec))/new_l1);
 	  new_aniso2 = 1 + par.lambda_aniso*(1 - fabs(InnerProduct((new_p-i_plus_1), c.cellvec))/new_l2);
 	  length_dh += w1 * ( new_aniso1 * DSQR(new_l1 - Node::target_length) -
 	  old_aniso1 * DSQR(old_l1 - Node::target_length) )
 	  +  w2 * ( new_aniso2 * DSQR(new_l2 - Node::target_length) -
 	  old_aniso2 * DSQR(old_l2 - Node::target_length) );
 	  } else {
 	  */
 	  length_dh += 2*Node::target_length * ( w1*(old_l1 - new_l1) +
 						 w2*(old_l2 - new_l2) ) +
 	    w1*(DSQR(new_l1)
 		- DSQR(old_l1))
 	    + w2*(DSQR(new_l2)
 		  - DSQR(old_l2));
 	  //}
+	}
 	// bending energy also holds for outer boundary
 	// first implementation. Can probably be done more efficiently
 	// calculate circumcenter radius (gives local curvature)
 	// the ideal bending state is flat... (K=0)
 	// if (cit->cell==-1 && node.cells.size()>2 /* boundary_pol, cell and a second cell */)  {
+	{
 	  // strong bending energy to resist "cleaving" by division planes
 	  double r1, r2, xc, yc;
 	  CircumCircle(i_min_1.x, i_min_1.y, old_p.x, old_p.y, i_plus_1.x, i_plus_1.y,
 		       &xc,&yc,&r1);
 	  CircumCircle(i_min_1.x, i_min_1.y, new_p.x, new_p.y, i_plus_1.x, i_plus_1.y,
 		       &xc,&yc, &r2);
 	  if (r1<0 || r2<0) {
 	    MyWarning::warning("r1 = %f, r2 = %f",r1,r2);
+	  }
 	  bending_dh += DSQR(1/r2 - 1/r1);
 	  //bending_dh += ( DSQR(1/r2) - DSQR(1/r1) );
 	  //cerr << "bending_dh = " << par.bend_lambda*bending_dh << endl;
+	}
 	/*cerr << "Bending = " <<  ( DSQR(1/r2) - DSQR(1/r1))  << endl;
 	  cerr << node.index << ": " << bending_dh << " (" << r1 << ", " << r2 << ") " << cit->nb1 << ", " << cit->nb2 << ")" << endl;
 	  }*/
 	/*else
 	  bending_dh += ( DSQR(1/r2) - DSQR(1/r1) );*/
 	/*if (cit->cell==-1) {
 	  cerr << node.index << ": " << bending_dh << " (" << r1 << ", " << r2 << ") " << cit->nb1 << ", " << cit->nb2 << ")" << endl;
 	  }*/
 	//bending_dh += r1 - r2;
+      }
       //const double bend_lambda = 100000;
       //const double bend_lambda = 0;
       /* double dh = //(area_dev_sum_after - area_dev_sum_before) +
 	 area_dh + par.lambda_celllength * cell_length_dh +
 	 par.lambda_length * length_dh + par.bend_lambda * bending_dh + par.alignment_lambda * alignment_dh;*/
       dh = //(area_dev_sum_after - area_dev_sum_before) +
 	area_dh + cell_length_dh +
 	par.lambda_length * length_dh + par.bend_lambda * bending_dh + par.alignment_lambda * alignment_dh;
       //cerr << "cell_length_dh = " << par.lambda_celllength * cell_length_dh << endl;
       //(length_constraint_after - length_constraint_before);
       if (node.fixed) {
 	// search the fixed cell to which this node belongs
 	// and displace these cells as a whole
 	// WARNING: undefined things will happen for connected fixed cells...
 	for (list<Neighbor>::iterator c=node.owners.begin();
 	     c!=node.owners.end();
 	     c++) {
 	for (list<Neighbor>::iterator c=node.owners.begin(); c!=node.owners.end(); c++) {
 	  if (!c->cell->BoundaryPolP() && c->cell->FixedP()) {
 	    sum_dh+=c->cell->Displace(rx,ry,0);
+	  }
+	}
       } else {
 	if (dh<-sum_stiff || RANDOM()<exp((-dh-sum_stiff)/par.T)) {
 	  // update areas of cells
 	  list<DeltaIntgrl>::const_iterator di_it = delta_intgrl_list.begin();
 	  for (list<Neighbor>::iterator cit=node.owners.begin();
 	       cit!=node.owners.end();
 	       ( cit++) ) {
 	    //m->getCell(cit->cell).area -= *da_it;
 	    //if (cit->cell>=0) {
 	  for (list<Neighbor>::iterator cit=node.owners.begin(); cit!=node.owners.end(); ( cit++) ) {
 	    if (!cit->cell->BoundaryPolP()) {
 	      cit->cell->area -= di_it->area;
 	      if (par.lambda_celllength) {
 		cit->cell->intgrl_x -= di_it->ix;
 		cit->cell->intgrl_y -= di_it->iy;
 		cit->cell->intgrl_xx -= di_it->ixx;
 		cit->cell->intgrl_xy -= di_it->ixy;
 		cit->cell->intgrl_yy -= di_it->iyy;
+	      }
 	      di_it++;
+	    }
+	  }
 	  double old_nodex, old_nodey;
 	  old_nodex=node.x;
 	  old_nodey=node.y;
 	  node.x = new_p.x;
 	  node.y = new_p.y;
 	  for (list<Neighbor>::iterator cit=node.owners.begin();
 	       cit!=node.owners.end();
 	       ( cit++) ) {
 	    /*   if (cit->cell >= 0 && cells[cit->cell].SelfIntersect()) {
 		 node.x = old_nodex;
 		 node.y = old_nodey;
 		 goto next_node;
 		 }*/
+	  }
 	  sum_dh += dh;
+	}
+      }
+    }
   next_node:
     delta_intgrl_list.clear();//dA_list.clear();
+  }
   return sum_dh;
+}
 void Mesh::InsertNode(Edge &e) {
   // Construct a new node in the middle of the edge
   Node *new_node = AddNode( new Node ( ( *e.first + *e.second )/2 ) );
   // if new node is inserted into the boundary
   // it will be part of the boundary, fixed, and source, too
   // The new node is part of the boundary only if both its neighbors are boundary nodes and the boundray proceeds from first to second.
   new_node->boundary = (e.first->BoundaryP() && e.first->BoundaryP()) && ((findNextBoundaryNode(e.first))->Index() == e.second->Index());
   new_node->fixed = e.first->fixed && e.second->fixed;
   new_node->sam = new_node->boundary && (e.first->sam || e.second->sam);
   // insert it into the boundary polygon;
   /* if (new_node->boundary) {
   // find the position of the first node in the boundary
   list<Node *>::iterator ins_pos = find
   (boundary_polygon->nodes.begin(),
   boundary_polygon->nodes.end(),
   e.first);
   // ... second node comes before or after it ...
   if (*(++ins_pos!=boundary_polygon->nodes.end()?
   ins_pos:boundary_polygon->nodes.begin())!=e.second) {
   boundary_polygon->nodes.insert(((ins_pos--)!=boundary_polygon->nodes.begin()?ins_pos:(--boundary_polygon->nodes.end())), new_node);
   // .. set the neighbors of the new node ...
   // in this case e.second and e.first are inverted
   new_node->owners.push_back( Neighbor(boundary_polygon, e.second, e.first ) );
   //cerr << "pushing back " << Neighbor(boundary_polygon->index, e.second, e.first ) << endl;
   } else {
   // insert before second node, so leave ins_pos as it is,
   // that is incremented
   boundary_polygon->nodes.insert(ins_pos, new_node);
   // .. set the neighbors of the new node ...
   new_node->owners.push_back( Neighbor(boundary_polygon, e.first, e.second ) );
   // cerr << "pushing back " << Neighbor(boundary_polygon->index, e.second, e.first ) << endl;
+  }
   }*/
   list<Neighbor> owners;
   // push all cells owning the two nodes of the divided edges
   // onto a list
   copy(e.first->owners.begin(),
        e.first->owners.end(),
        back_inserter(owners));
@@ @@ -1142,1128 +1054,977 @@ void Mesh::InsertNode(Edge &e) { @@
       if (cpos->nb1 == e.second) {
 	cpos->nb1 = new_node;
       } else
 	if (cpos->nb2 == e.second) {
 	  cpos->nb2 = new_node;
+	}
       // - same for Node e.second
       cpos=
 	find_if( e.second->owners.begin(),
 		 e.second->owners.end(),
 		 bind2nd( mem_fun_ref(&Neighbor::CellEquals), c->cell->Index()) );
       // - correct the record
       if (cpos->nb1 == e.first) {
 	cpos->nb1 = new_node;
       } else
 	if (cpos->nb2 == e.first) {
 	  cpos->nb2 = new_node;
+	}
     } else break;
     c++;
+  }
   // Repair neighborhood lists in a second loop, to make sure all
   // `administration' is up to date
   while (c!=owners.end()) {
     c=adjacent_find(c,owners.end(),  neighbor_cell_eq);
     // repair neighborhood lists of cell and Wall lists
     //if (c->cell>=0) {
     if (!c->cell->BoundaryPolP()) {
       c->cell->ConstructNeighborList();
       //      debug_stream << "Repairing NeighborList of " << c->cell << endl;
+    }
     c++;
+  }
   //  debug_stream.flush();
+}
 /*
   Calculate circumcircle of triangle (x1,y1), (x2,y2), (x3,y3)
   The circumcircle centre is returned in (xc,yc) and the radius in r
   NOTE: A point on the edge is inside the circumcircle
 */
 void Mesh::CircumCircle(double x1,double y1,double x2,double y2,double x3,double y3,
 			double *xc,double *yc,double *r)
+{
   double m1,m2,mx1,mx2,my1,my2;
   double dx,dy,rsqr;
   /* Check for coincident points */
   /*if (abs(y1-y2) < TINY && abs(y2-y3) < TINY)
     return(false);*/
   if (abs(y2-y1) < TINY) {
     m2 = - (x3-x2) / (y3-y2);
     mx2 = (x2 + x3) / 2.0;
     my2 = (y2 + y3) / 2.0;
     *xc = (x2 + x1) / 2.0;
     *yc = m2 * (*xc - mx2) + my2;
   } else if (abs(y3-y2) < TINY) {
     m1 = - (x2-x1) / (y2-y1);
     mx1 = (x1 + x2) / 2.0;
     my1 = (y1 + y2) / 2.0;
     *xc = (x3 + x2) / 2.0;
     *yc = m1 * (*xc - mx1) + my1;
   } else {
     m1 = - (x2-x1) / (y2-y1);
     m2 = - (x3-x2) / (y3-y2);
     mx1 = (x1 + x2) / 2.0;
     mx2 = (x2 + x3) / 2.0;
     my1 = (y1 + y2) / 2.0;
     my2 = (y2 + y3) / 2.0;
     *xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
     *yc = m1 * (*xc - mx1) + my1;
+  }
   dx = x2 - *xc;
   dy = y2 - *yc;
   rsqr = dx*dx + dy*dy;
   *r = sqrt(rsqr);
   return;
   // Suggested
   // return((drsqr <= rsqr + EPSILON) ? TRUE : FALSE);
+}
 //
 // return the total amount of chemical "ch" in the leaf
 double Mesh::SumChemical(int ch) {
   double sum=0.;
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
   for (vector<Cell *>::iterator i=cells.begin(); i!=cells.end(); i++) {
     sum+=(*i)->chem[ch];
+  }
   return sum;
+}
 void Mesh::CleanUpCellNodeLists(void) {
   typedef vector <vector<Cell *>::iterator> CellItVect;
   CellItVect cellstoberemoved;
   vector<int> cellind;
   // Start of by removing all stale walls.
   //DeleteLooseWalls();
   // collect all dead cells that need to be removed from the simulation
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
   for (vector<Cell *>::iterator i=cells.begin(); i!=cells.end(); i++) {
     if ((*i)->DeadP()) {
       // collect the iterators
       cellstoberemoved.push_back(i);
       // collect the indices
       cellind.push_back((*i)->index);
     } else {
       // Remove pointers to dead Walls
       for (list<Wall *>::iterator w=(*i)->walls.begin();
 	   w!=(*i)->walls.end();
 	   w++) {
       for (list<Wall *>::iterator w=(*i)->walls.begin(); w!=(*i)->walls.end(); w++) {
 	if ((*w)->DeadP()) {
 	  (*w)=0;
+	}
+      }
       (*i)->walls.remove(0);
+    }
+  }
   // Remove pointers to dead Walls from BoundaryPolygon
   for (list<Wall *>::iterator w=boundary_polygon->walls.begin();
        w!=boundary_polygon->walls.end();
        w++) {
   for (list<Wall *>::iterator w=boundary_polygon->walls.begin(); w!=boundary_polygon->walls.end(); w++) {
     if ((*w)->DeadP()) {
       (*w)=0;
+    }
+  }
   boundary_polygon->walls.remove(0);
   // Renumber cells; this is most efficient if the list of dead cell indices is sorted
   sort(cellind.begin(),cellind.end());
   // Reindexing of Cells
   for (vector<int>::reverse_iterator j=cellind.rbegin();
        j!=cellind.rend();
        j++) {
     for (vector<Cell *>::reverse_iterator i=cells.rbegin();
 	 i!=cells.rend();
 	 i++) {
   for (vector<int>::reverse_iterator j=cellind.rbegin(); j!=cellind.rend(); j++) {
     for (vector<Cell *>::reverse_iterator i=cells.rbegin(); i!=cells.rend(); i++) {
       if (*j < (*i)->index) (*i)->index--;
+    }
+  }
   // Actual deleting of Cells
   // We must delete in reverse order, otherwise the iterators become redefined
   for ( CellItVect::reverse_iterator i=cellstoberemoved.rbegin();
 	i!=cellstoberemoved.rend();
 	i++) {
   for ( CellItVect::reverse_iterator i=cellstoberemoved.rbegin(); i!=cellstoberemoved.rend(); i++) {
     Cell::NCells()--;
     cells.erase(*i);
+  }
   // same for nodes
   typedef vector <vector<Node *>::iterator> NodeItVect;
   NodeItVect nodestoberemoved;
   vector<int> nodeindlist;
   // collect iterators and indices of dead nodes
   for (vector<Node *>::iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (vector<Node *>::iterator i=nodes.begin(); i!=nodes.end(); i++) {
     if ((*i)->DeadP()) {
       nodestoberemoved.push_back( i );
       nodeindlist.push_back((*i)->index);
+    }
+  }
   // sort the list of dead nodes for renumbering
   sort(nodeindlist.begin(),nodeindlist.end());
   // Reindicing of Nodes
   for (vector<int>::reverse_iterator j=nodeindlist.rbegin();
        j!=nodeindlist.rend();
        j++) {
     for (vector<Node *>::reverse_iterator i=nodes.rbegin();
 	 i!=nodes.rend();
 	 i++) {
   for (vector<int>::reverse_iterator j=nodeindlist.rbegin(); j!=nodeindlist.rend(); j++) {
     for (vector<Node *>::reverse_iterator i=nodes.rbegin(); i!=nodes.rend(); i++) {
       if (*j < (*i)->index) {
 	(*i)->index--;
+      }
+    }
+  }
   // Actual deleting of nodes
   // We must delete in reverse order, otherwise the iterators become redefined
   for ( NodeItVect::reverse_iterator i=nodestoberemoved.rbegin();
 	i!=nodestoberemoved.rend();
 	i++) {
   for ( NodeItVect::reverse_iterator i=nodestoberemoved.rbegin(); i!=nodestoberemoved.rend(); i++) {
     Node::nnodes--;
     nodes.erase(*i);
+  }
   for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
   for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
     if ((*w)->DeadP()) {
       Wall::nwalls--;
       delete *w;
       *w = 0;
+    }
+  }
   walls.remove( 0 );
   // Clean up all intercellular connections and redo everything
   for (vector<Node *>::iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (vector<Node *>::iterator i=nodes.begin(); i!=nodes.end(); i++) {
     (*i)->owners.clear();
+  }
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
   for (vector<Cell *>::iterator i=cells.begin(); i!=cells.end(); i++) {
     (*i)->ConstructConnections();
+  }
   boundary_polygon->ConstructConnections();
   /* for (list<Wall *>::iterator w=walls.begin();
      w!=walls.end();
      w++) {
      delete *w;
+     }
      walls.clear();
      cerr << "Cleared walls\n";
      for (vector<Cell *>::iterator i=cells.begin();
      i!=cells.end();
      i++) {
      (*i)->ConstructWalls();
+     }
   */
   // remake shuffled_nodes and shuffled cells
   shuffled_nodes.clear();
   shuffled_nodes = nodes;
   shuffled_cells.clear();
   shuffled_cells = cells;
+}
 void Mesh::CutAwayBelowLine( Vector startpoint, Vector endpoint) {
   // Kills all cells below the line startpoint -> endpoint
   Vector perp = (endpoint-startpoint).Perp2D().Normalised();
 #ifdef QDEBUG
   qDebug() << "Before Apoptose" << endl;
 #endif
   TestIllegalWalls();
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
   for (vector<Cell *>::iterator i=cells.begin(); i!=cells.end(); i++) {
     // do some vector geometry to check whether the cell is below the cutting line
     Vector cellvec = ((*i)->Centroid()-startpoint);
     if ( InnerProduct(perp, cellvec) < 0 ) {
       // remove those cells
       (*i)->Apoptose();
+    }
+  }
 #ifdef QDEBUG
   qDebug() << "Before CleanUpCellNodeLists" << endl;
 #endif
   TestIllegalWalls();
   CleanUpCellNodeLists();
+}
 void Mesh::CutAwaySAM(void) {
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
   for (vector<Cell *>::iterator i=cells.begin(); i!=cells.end(); i++) {
     if( (*i)->Boundary() == Cell::SAM ) {
       (*i)->Apoptose();
+    }
+  }
   TestIllegalWalls();
   CleanUpCellNodeLists();
+}
 void Mesh::TestIllegalWalls(void) {
   for (list<Wall *>::iterator w = walls.begin();
        w!=walls.end();
        w++) {
   for (list<Wall *>::iterator w = walls.begin(); w!=walls.end(); w++) {
     if ((*w)->IllegalP() ) {
 #ifdef QDEBUG
       qDebug() << "Wall " << **w << " is illegal." << endl;
 #endif
+    }
+  }
+}
 class node_owners_eq : public unary_function<Node, bool> {
   int no;
 public:
   explicit node_owners_eq(int nn) { no=nn; }
   bool operator() (const Node &n) const {
     if (n.CellsSize()==1)
       return true;
     else
       return false;
+  }
 };
 void Mesh::RepairBoundaryPolygon(void) {
   // After serious manipulations (e.g. after cutting part off the
   // leaf) repair the boundary polygon. It assumes the cut line has
   // already been marked "boundary" and the original boundary marks
   // were not removed.
   //
   // So, this function just puts boundary nodes into the boundary
   // polygon in the right order; it cannot detect boundaries from
   // scratch.
   Node *boundary_node=0, *next_boundary_node=0, *internal_node;
   set<int> original_boundary_nodes, repaired_boundary_nodes;
   vector<int> difference; // set difference result
   // Step 0: print out current boundary polygon
 #ifdef QDEBUG
   qDebug() << endl << "Original Boundary Polygon node indices: ";
   foreach (Node* node, boundary_polygon->nodes) {
     qDebug() << node->Index() << " " ;
+  }
   qDebug() << endl << endl;
 #endif
   // Step 1a: Create a set containing the current boundary polygon nodes' Indices.
   foreach (Node* node, boundary_polygon->nodes) {
     original_boundary_nodes.insert(node->Index());
+  }
   // Step 1b: remove all nodes from boundary polygon
   boundary_polygon->nodes.clear();
   // Step 2: Remove all references to the boundary polygon from the Mesh's current list of nodes
   foreach (Node* node, nodes) {
     node->Unmark(); // remove marks, we need them to determine if we have closed the circle
     list<Neighbor>::iterator boundary_ref_pos;
     if ((boundary_ref_pos = find_if (node->owners.begin(), node->owners.end(),
 				     bind2nd(mem_fun_ref(&Neighbor::CellEquals), -1))) != node->owners.end()) {
       // i.e. if one of the node's owners is the boundary polygon
       node->owners.erase(boundary_ref_pos); // remove the reference
+    }
+  }
   // Step 3: Search for the first boundary node.  We reconstruct the
   // boundary polygon by moving along the boundary nodes until we've
   // encircled the polygon. Since manually adding nodes may have
   // turned nodes previously along the boundary into internal nodes,
   // we search through all the node until we find first boundary node
   // and proceed from there. If findNextBoundaryNode() returns a node
   // other than the one passed to it, the original node is the first
   // boundary node.
   foreach (Node* node, nodes) {
     if ((findNextBoundaryNode(node))->index != node->index){
       next_boundary_node = node;
       break;
+    }
+  }
   // We have a problem if we arrive here without having found a boundary node.
   if (!next_boundary_node) throw("Cannot find a boundary node!.");
   // Reconstruct the list of boundary polygon nodes.
   do {
     boundary_node = next_boundary_node;
     boundary_node->Mark();
     boundary_polygon->nodes.push_back(boundary_node);
     next_boundary_node = findNextBoundaryNode(boundary_node);
   } while ( !next_boundary_node->Marked() );
   // Create a set containing the reconstructed boundary polygon nodes' Indices.
   for (list<Node *>::iterator it = boundary_polygon->nodes.begin(); it!=boundary_polygon->nodes.end(); ++it) {
     repaired_boundary_nodes.insert((*it)->Index());
+  }
   // Calculate the difference between the original and repaired sets of boundary nodes
   // yielding the set of nodes that are no longer part of the boundary polygon.
   set_difference(original_boundary_nodes.begin(), original_boundary_nodes.end(),
                  repaired_boundary_nodes.begin(), repaired_boundary_nodes.end(), back_inserter(difference));
   // Tell each node in the difference that it's no longer part of the boundary polygon
   vector<Node *>::iterator internal_node_it;
   foreach (int i, difference){
     internal_node_it = find_if (nodes.begin(), nodes.end(), bind2nd(mem_fun(&Node::IndexEquals), i));
     internal_node = *internal_node_it; // dereference the itterator to get to the node pointer
     if (!internal_node) throw("Found a null Node pointer.");
     internal_node->UnsetBoundary();
+  }
   boundary_polygon->ConstructConnections();
   for (list<Wall *>::iterator w=boundary_polygon->walls.begin();
        w!=boundary_polygon->walls.end();
        w++) {
   for (list<Wall *>::iterator w=boundary_polygon->walls.begin(); w!=boundary_polygon->walls.end(); w++) {
     if ((*w)->DeadP()) {
       (*w)=0;
+    }
+  }
   boundary_polygon->walls.remove(0);
   boundary_polygon->ConstructNeighborList();
 #ifdef QDEBUG
   qDebug() << "Repaired Boundary Polygon node indices: ";
   foreach (Node* node, boundary_polygon->nodes){
     qDebug() << node->Index() << " " ;
+  }
   qDebug() << endl ;
 #ifdef _undefined_
   qDebug() << "NODES:" << endl;
   foreach(Node* node, nodes) {
     qDebug() << *node;
+  }
   qDebug() << endl;
   qDebug() << "WALLS:" << endl;
   foreach(Wall* wall, walls) {
     qDebug() << *wall;
+  }
   qDebug() << endl;
   qDebug() << "CELLS:" << endl;
   foreach(Cell* cell, cells) {
     qDebug() << *cell;
+  }
   qDebug() << endl;
 #endif
 #endif
+}
 Node* Mesh::findNextBoundaryNode(Node* boundary_node) {
   bool found_next_boundary_node = false;
   Node *next_boundary_node = NULL;
   set<int> boundary_node_owners; // This is a list of the current boundary node's owners' Ids
   vector<int> neighborIds; // A list of the current boundary node's owners' 2nd neighbor Ids
   vector<set<int> *>  nodeOwners; // A vector of set pointers where each set contains the owner Ids of the nodes in the neighborIds list.
   vector<int> intersection; // set intersection result
   // The next boundary node is that which has only one owner in common with the current boundary node
   for (list<Neighbor>::iterator it=boundary_node->owners.begin(); it!=boundary_node->owners.end(); ++it) {
     if (it->cell->Index() != -1) boundary_node_owners.insert(it->cell->Index()); // Save each of the current boundary node's owners' Ids - except the boundary polygon
     set<int> *owners = new set<int>; // create a set to hold a 2nd neighbor's owners' Ids
     nodeOwners.push_back(owners);
     neighborIds.push_back(it->nb2->Index());
     foreach(Neighbor neighbor, it->nb2->owners){
       if (neighbor.cell->Index() != -1) owners->insert(neighbor.cell->Index()); // Save second neighbors' owners' Ids - except the boundary polygon
+    }
+  }
   vector<int>::iterator itt = neighborIds.begin();
   vector<set<int> *>::iterator it = nodeOwners.begin();
 #ifdef QDEBUG
   qDebug() << "Boundary node: " <<  boundary_node->Index() << " is owned by the following cells: ";
   foreach (int i, boundary_node_owners){
     qDebug() << i << "  ";
+  }
   qDebug() << endl;
 #endif
   for (; it < nodeOwners.end(); it++, itt++) {
     intersection.clear();
     set_intersection(boundary_node_owners.begin(), boundary_node_owners.end(), (*it)->begin(), (*it)->end(), back_inserter(intersection));
 #ifdef QDEBUG
     qDebug() << "The intersection of the boundary node(" << boundary_node->Index() << ") owners and its 2nd neighbor(" <<  *itt << ") owners is: ";
     foreach (int i, intersection){
       qDebug() << i << "  ";
+    }
     qDebug() << endl;
 #endif
     if (intersection.size() == 1){
       found_next_boundary_node = true;
       vector<Node *>::iterator next_boundary_node_it = find_if (nodes.begin(), nodes.end(), bind2nd(mem_fun(&Node::IndexEquals), *itt));
       next_boundary_node = *next_boundary_node_it; // defeference the itterator to get to the node pointer
 #ifdef QDEBUG
       qDebug() << "The Current boundary node is: " << boundary_node->Index()
 	       << ". The Next boundary node is: " << *itt << ((next_boundary_node->Marked()) ? " Marked" : " Unmarked") << endl << endl;
 #endif
       break;
+    }
+  }
 #ifdef QDEBUG
   if (!found_next_boundary_node) {
     qDebug() << "OOPS! Didn't find the next boundrary node!" << endl;
+  }
 #endif
   return next_boundary_node;
+}
 void Mesh::CleanUpWalls(void) {
   for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
   for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
     if ((*w)->DeadP()) {
       delete *w;
       (*w)=0;
+    }
+  }
   walls.remove(0);
+}
 void Mesh::Rotate(double angle, Vector center) {
   // Rotate the mesh over the angle "angle", relative to center point "center".
   Matrix rotmat;
   rotmat.Rot2D(angle);
   for (vector<Node *>::iterator n=nodes.begin();
        n!=nodes.end();
        n++) {
   for (vector<Node *>::iterator n=nodes.begin(); n!=nodes.end(); n++) {
     (*n)->setPos ( rotmat * ( *(*n) - center ) + center );
+  }
+}
 void Mesh::PrintWallList( void ) {
   transform ( walls.begin(), walls.end(), ostream_iterator<Wall>(cerr, "\n"), deref_ptr<Wall> );
+}
 #include <QString>
 //#include "forwardeuler.h"
 #include "rungekutta.h"
 class SolveMesh : public RungeKutta {
 private:
   SolveMesh(void);
 public:
   SolveMesh(Mesh *m_) {
     m = m_;
     kmax=0;
     kount=0;
     xp=0; yp=0; dxsav=0;
+  }
 protected:
   virtual void derivs(double x, double *y, double *dydx) {
     // set mesh with new values given by ODESolver
     // (we must do this, because only mesh knows the connections
     // between the variables)
     m->setValues(x,y);
     m->Derivatives(dydx);
     /*static int c=0;
       QString fname("derivs%1.dat");
       ofstream of(fname.arg(++c).ascii());
       for (int i=0;i<m->NEqs();i++) {
       of << x << " " << dxdy[i] << endl;
+      }
       of.close();
     */
     //cerr << "Calculated derivatives at " << x << "\n";
+  }
 private:
   Mesh *m;
   int kmax,kount;
   double *xp,**yp,dxsav;
   bool monitor_window;
 };
 void Mesh::ReactDiffuse(double delta_t) {
   // Set Lengths of Walls
   for_each ( walls.begin(), walls.end(),
 	     mem_fun( &Wall::SetLength ) );
   static SolveMesh *solver = new SolveMesh(this);
   int nok, nbad, nvar;
   double *ystart = getValues(&nvar);
   solver->odeint(ystart, nvar, getTime(), getTime() + delta_t,
 		 par.ode_accuracy, par.dt, 1e-10, &nok, &nbad);
   setTime(getTime()+delta_t);
   setValues(getTime(),ystart);
+}
 Vector Mesh::FirstConcMoment(int chem) {
   Vector moment;
   for (vector<Cell *>::const_iterator c=cells.begin();
        c!=cells.end();
        c++) {
   for (vector<Cell *>::const_iterator c=cells.begin(); c!=cells.end(); c++) {
     moment += (*c)->Chemical(chem) * (*c)->Centroid();
+  }
   return moment / (double)cells.size();
+}
 /*! This member function deletes all walls connected to two dead cells from the mesh.
   It should be called before the Cells are actually removed.
   If the cell is connect to one dead cell only, that reference is substituted for a reference
   to the boundary polygon.
 */
 void Mesh::DeleteLooseWalls(void) {
   list<Wall *>::iterator w=walls.begin();
   while (w!=walls.end()) {
     // if both cells of the wall are dead, remove the wall
     if ((*w)->C1()->DeadP() || (*w)->C2()->DeadP()) {
       if ((*w)->C1()->DeadP() && (*w)->C2()->DeadP()) {
 	delete *w;
 	w=walls.erase(w);
       } else {
 	if ((*w)->C1()->DeadP())
 	  (*w)->c1 = boundary_polygon;
 	else
 	  (*w)->c2 = boundary_polygon;
 	w++;
+      }
     } else {
       w++;
+    }
+  }
+}
 /*void Mesh::FitLeafToCanvas(double width, double height) {
   Vector bbll,bbur;
   BoundingBox(bbll,bbur);
   double scale_x = width/(bbur.x-bbll.x);
   double scale_y = height/(bbur.y-bbll.y);
   double factor = scale_x<scale_y ? scale_x:scale_y;
   Cell::SetMagnification(factor); // smallest of scale_x and scale_y
   double offset_x = (width/Cell::Magnification()-(bbur.x-bbll.x))/2.;
   double offset_y = (height/Cell::Magnification()-(bbur.y-bbll.y))/2.;
   Cell::setOffset(offset_x, offset_y);
   }*/
 void Mesh::CleanChemicals(const vector<double> &clean_chem) {
   if (clean_chem.size()!=(unsigned)Cell::NChem()) {
     throw "Run time error in Mesh::CleanChemicals: size of clean_chem should be equal to Cell::NChem()";
+  }
   for (vector<Cell *>::iterator c=cells.begin();
        c!=cells.end();
        c++) {
   for (vector<Cell *>::iterator c=cells.begin(); c!=cells.end(); c++) {
     for (int i=0;i<Cell::NChem();i++) {
       (*c)->SetChemical(i,clean_chem[i]);
+    }
     (*c)->SetNewChemToChem();
+  }
+}
 void Mesh::CleanTransporters(const vector<double> &clean_transporters) {
   if (clean_transporters.size()!=(unsigned)Cell::NChem()) {
     throw "Run time error in Mesh::CleanTransporters: size ofclean_transporters should be equal to Cell::NChem()";
+  }
   // clean transporters
   for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
   for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
     for (int i=0;i<Cell::NChem();i++) {
       (*w)->setTransporters1(i,clean_transporters[i]); (*w)->setNewTransporters1(i,clean_transporters[i]);
       (*w)->setTransporters2(i,clean_transporters[i]); (*w)->setNewTransporters2(i,clean_transporters[i]);
+    }
+  }
+}
 void Mesh::RandomizeChemicals(const vector<double> &max_chem, const vector<double> &max_transporters) {
   if (max_chem.size()!=(unsigned)Cell::NChem() || max_transporters.size()!=(unsigned)Cell::NChem()) {
     throw "Run time error in Mesh::CleanChemicals: size of max_chem and max_transporters should be equal to Cell::NChem()";
+  }
   for (vector<Cell *>::iterator c=cells.begin();
        c!=cells.end();
        c++) {
   for (vector<Cell *>::iterator c=cells.begin(); c!=cells.end(); c++) {
     for (int i=0;i<Cell::NChem();i++) {
       (*c)->SetChemical(i,max_chem[i]*RANDOM());
+    }
     (*c)->SetNewChemToChem();
+  }
   // randomize transporters
   for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
   for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
     for (int i=0;i<Cell::NChem();i++) {
       (*w)->setTransporters1(i,max_transporters[i] * RANDOM()); (*w)->setNewTransporters1(i, (*w)->Transporters1(i) );
       (*w)->setTransporters2(i,max_transporters[i] * RANDOM()); (*w)->setNewTransporters2(i, (*w)->Transporters1(i) );
+    }
+  }
+}
 //!\brief Calculates a vector with derivatives of all variables, which
 // we can pass to an ODESolver.
 void Mesh::Derivatives(double *derivs) {
   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;
   //static double *derivs = 0;
   // derivs is allocated by RungeKutta class.
   for (int i=0;i<neqs;i++) { derivs[i]=0.;}
   for (int i=0;i<neqs;i++) {
     derivs[i]=0.;
+  }
   // Layout of derivatives: cells [ chem1 ... chem n]  walls [ [ w1(chem 1) ... w1(chem n) ] [ w2(chem 1) ... w2(chem n) ] ]
   int i=0;
   for (vector<Cell *>::iterator c=cells.begin();
        c!=cells.end();
        c++) {
     //(*cr)(*c, &(derivs[i]));
   for (vector<Cell *>::iterator c=cells.begin(); c!=cells.end(); c++) {
     plugin->CellDynamics(*c, &(derivs[i]));
     i+=nchems;
+  }
   for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
   for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
     // (*wr)(*w, &(derivs[i]), &(derivs[i+nchems]));
     plugin->WallDynamics(*w,  &(derivs[i]), &(derivs[i+nchems]));
     // Transport function adds to derivatives of cell chemicals
     double *dchem_c1 = &(derivs[(*w)->c1->Index() * nchems]);
     double *dchem_c2 = &(derivs[(*w)->c2->Index() * nchems]);
     //plugin->CelltoCellTransport(*w, &(derivs[(*w)->c1->Index() * nchems]),
     //	  &(derivs[(*w)->c2->Index() * nchems]));
     // quick fix: dummy values to prevent end user from writing into outer space and causing a crash :-)
     // start here if you want to implement chemical input/output into environment over boundaries
     double dummy1, dummy2;
     if ((*w)->c1->Index()<0) { // tests if c1 is the boundary pol
       dchem_c1 = &dummy1;
+    }
     if ((*w)->c2->Index()<0) {
       dchem_c2 = &dummy2;
+    }
     plugin->CelltoCellTransport(*w, dchem_c1, dchem_c2);
     //(*tf)(*w, &(derivs[(*w)->c1->Index() * nchems]),
     //&(derivs[(*w)->c2->Index() * nchems] ) );
     i+=2*nchems;
+  }
+}
 void Mesh::setValues(double x, double *y) {
   //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;
   // Layout of derivatives: cells [ chem1 ... chem n]  walls [ [ w1(chem 1) ... w1(chem n) ] [ w2(chem 1) ... w2(chem n) ] ]
   int i=0;
   static int emit_count=0;
   const int stride = 100;
   for (vector<Cell *>::iterator c=cells.begin();
        c!=cells.end();
        c++) {
   for (vector<Cell *>::iterator c=cells.begin(); c!=cells.end(); c++) {
     for (int ch=0;ch<nchems;ch++) {
       (*c)->SetChemical(ch, y[i+ch]);
+    }
     if ( !(emit_count%stride)) {
       (*c)->EmitValues(x);
+    }
     i+=nchems;
+  }
   for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
   for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
     for (int ch=0;ch<nchems;ch++) {
       (*w)->setTransporters1(ch,y[i+ch]);
+    }
     i+=nchems;
     for (int ch=0;ch<nchems;ch++) {
       (*w)->setTransporters2(ch,y[i+ch]);
+    }
     i+=nchems;
+  }
   emit_count++;
+}
 double *Mesh::getValues(int *neqs) {
   int nwalls = walls.size();
   int ncells = cells.size();
   int nchems = Cell::NChem();
   // two eqs per chemical for each wall, 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.
   (*neqs) = 2 * nwalls * nchems + ncells * nchems;
   // Layout of derivatives: cells [ chem1 ... chem n]  walls [ [ w1(chem 1) ... w1(chem n) ] [ w2(chem 1) ... w2(chem n) ] ]
   static double *values = 0;
   if (values!=0) { delete[] values; }
   values = new double[*neqs];
   int i=0;
   for (vector<Cell *>::iterator c=cells.begin();
        c!=cells.end();
        c++) {
   for (vector<Cell *>::iterator c=cells.begin(); c!=cells.end(); c++) {
     for (int ch=0;ch<nchems;ch++) {
       values[i+ch]=(*c)->Chemical(ch);
+    }
     i+=nchems;
+  }
   for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
   for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
     for (int ch=0;ch<nchems;ch++) {
       values[i+ch]=(*w)->Transporters1(ch);
+    }
     i+=nchems;
     for (int ch=0;ch<nchems;ch++) {
       values[i+ch]=(*w)->Transporters2(ch);
+    }
     i+=nchems;
+  }
   return values;
+}
 void Mesh::DrawNodes(QGraphicsScene *c) const {
   for (vector<Node *>::const_iterator n=nodes.begin();
        n!=nodes.end();
        n++) {
   for (vector<Node *>::const_iterator n=nodes.begin(); n!=nodes.end(); n++) {
     Node *i=*n;
     NodeItem *item = new NodeItem ( &(*i), c );
     item->setColor();
     item->setZValue(5);
     item->show();
     item ->setPos(((Cell::offset[0]+i->x)*Cell::factor),
 		  ((Cell::offset[1]+i->y)*Cell::factor) );
+  }
+}
 /*! Returns the sum of protein "ch" of a cycling protein in cells and walls */
 double Mesh::CalcProtCellsWalls(int ch) const {
   double sum_prot=0.;
   // At membranes
   for (list<Wall *>::const_iterator w=walls.begin();
        w!=walls.end();
        w++) {
   for (list<Wall *>::const_iterator w=walls.begin(); w!=walls.end(); w++) {
     sum_prot += (*w)->Transporters1(ch);
     sum_prot += (*w)->Transporters2(ch);
+  }
   // At cells
   for (vector<Cell *>::const_iterator c=cells.begin();
        c!=cells.end();
        c++) {
   for (vector<Cell *>::const_iterator c=cells.begin(); c!=cells.end(); c++) {
     sum_prot += (*c)->Chemical(ch);
+  }
   return sum_prot;
+}
 void Mesh::SettoInitVals(void) {
   vector<double> clean_chem(Cell::NChem());
   vector<double> clean_transporters(Cell::NChem());
   for (int i=0;i<Cell::NChem();i++) {
     clean_transporters[i]=0.;
     clean_chem[i]=par.initval[i];
+  }
   // Amount of PIN1
   //clean_chem[1] = 0.;
   CleanChemicals(clean_chem);
   CleanTransporters(clean_transporters);
+}
 string Mesh::getTimeHours(void) const {
   int hours = (int)(time / 3600);
   int mins = (int)((time - hours * 3600)/60);
   int secs = (int)((time - hours * 3600 - mins * 60));
   ostringstream tstr;
   tstr << hours << " h " << mins << " m " << secs << " s";
   return tstr.str();
+}
 QVector<qreal> Mesh::VertexAngles(void) {
   QVector<qreal> angles;
   for (vector<Node *>::const_iterator n=nodes.begin();
        n!=nodes.end();
        n++) {
   for (vector<Node *>::const_iterator n=nodes.begin(); n!=nodes.end(); n++) {
     if ((*n)->Value()>2 && !(*n)->BoundaryP() ) {
       angles+=(*n)->NeighbourAngles();
+    }
+  }
   return angles;
+}
 QVector< QPair<qreal,int> > Mesh::VertexAnglesValues(void) {
   QVector< QPair<qreal,int> > anglesvalues;
   for (vector<Node *>::const_iterator n=nodes.begin();
        n!=nodes.end();
        n++) {
   for (vector<Node *>::const_iterator n=nodes.begin(); n!=nodes.end(); n++) {
     if ((*n)->Value()>2 && !(*n)->BoundaryP() ) {
       QVector<qreal> angles = (*n)->NeighbourAngles();
       int value_vertex = angles.size();
       for (QVector<qreal>::ConstIterator i=angles.begin();
 	   i!=angles.end();
 	   i++) {
       for (QVector<qreal>::ConstIterator i=angles.begin(); i!=angles.end(); i++) {
 	anglesvalues += QPair< qreal, int > (*i, value_vertex);
+      }
+    }
+  }
   return anglesvalues;
+}
 void Mesh::Clean(void) {
 #ifdef QDEBUG
   qDebug() << "Freeing nodes" << endl;
 #endif
   for (vector<Node *>::iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (vector<Node *>::iterator i=nodes.begin(); i!=nodes.end(); i++) {
     delete *i;
+  }
   nodes.clear();
   Node::nnodes=0;
 #ifdef QDEBUG
   qDebug() << "Freeing node sets" << endl;
 #endif
   for (vector<NodeSet *>::iterator i=node_sets.begin();
        i!=node_sets.end();
        i++) {
   for (vector<NodeSet *>::iterator i=node_sets.begin(); i!=node_sets.end(); i++) {
     delete *i;
+  }
   node_sets.clear();
 #ifdef QDEBUG
   qDebug() << "Freeing cells" << endl;
 #endif
   //CellsStaticDatamembers *old_static_data_mem = Cell::GetStaticDataMemberPointer();
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
   for (vector<Cell *>::iterator i=cells.begin(); i!=cells.end(); i++) {
     delete *i;
+  }
   //Cell::static_data_members = old_static_data_mem;
   cells.clear();
   Cell::NCells()=0;
   delete boundary_polygon; // (already deleted during cleaning of cells?)
 #ifdef QDEBUG
   qDebug() << "Freeing walls" << endl;
 #endif
   for (list<Wall *>::iterator i=walls.begin();
        i!=walls.end();
        i++) {
   for (list<Wall *>::iterator i=walls.begin(); i!=walls.end(); i++) {
     delete *i;
+  }
   walls.clear();
   Wall::nwalls=0;
   node_insertion_queue.clear();
   shuffled_nodes.clear();
   shuffled_cells.clear();
   time = 0.0;
+}
 void Mesh::StandardInit(void) {
   boundary_polygon = new BoundaryPolygon();
   Cell &circle=CircularCell(0,0,10,10);
   circle.SetTargetArea(circle.CalcArea());
   circle.SetTargetLength(par.target_length);
   circle.SetLambdaLength(par.lambda_celllength);
   SetBaseArea();
   // clean up chemicals
   for (int c=0; c<Cell::NChem(); c++) {
     circle.SetChemical(c, 0.);
+  }
+}
+/* finis */

src/modelcatalogue.h

➞

Show inline comments

 /*
+ *
  *  $Id$
+ *
  *  This file is part of the Virtual Leaf.
+ *
  *  The Virtual Leaf is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
+ *
  *  The Virtual Leaf is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU General Public License for more details.
+ *
  *  You should have received a copy of the GNU General Public License
  *  along with the Virtual Leaf.  If not, see <http://www.gnu.org/licenses/>.
+ *
  *  Copyright 2010 Roeland Merks.
+ *
  */
 #ifndef _MODELCATALOGUE_H_
 #define _MODELCATALOGUE_H_
 #include "warning.h"
 #include "cell.h"
 #include "mesh.h"
 #include "simplugin.h"
 #include "canvas.h"
 #include <QPluginLoader>
 #include <Q3PopupMenu>
 #include <QAction>
 #include <QDir>
 #include <QApplication>
 #include <QObject>
 #include <QAction>
 #include <QMenu>
 class ModelCatalogue : public QObject {
   Q_OBJECT
     public:
   ModelCatalogue(Mesh *mesh, Main *mainwin, const char *model);
   void LoadPlugins();
   void LoadPlugin(const char *model);
   void InstallFirstModel();
   void PopulateModelMenu();
   public slots:
   void InstallModel(SimPluginInterface *model);
   void InstallModel(QAction *modelaction);
  private:
   QVector<SimPluginInterface *> models;
   Mesh *mesh;
   Main *mainwin;
 };
 #endif
 /* finis */

src/modelelement.h

➞

Show inline comments

 /*
+ *
  *  $Id$
+ *
  *  This file is part of the Virtual Leaf.
+ *
  *  The Virtual Leaf is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
+ *
  *  The Virtual Leaf is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU General Public License for more details.
+ *
  *  You should have received a copy of the GNU General Public License
  *  along with the Virtual Leaf.  If not, see <http://www.gnu.org/licenses/>.
+ *
  *  Copyright 2010 Roeland Merks.
+ *
  */
 #ifndef _MODELELEMENT_H_
 #define _MODELELEMENT_H_
 #include "vector.h"
 // pure virtual base class for all model elements
 class ModelElement {
  public:
   virtual void Displace(void) {};
   virtual ~ModelElement() {};
 };
 #endif
 /* finis */

src/node.cpp

➞

Show inline comments

@@ @@ -80,139 +80,134 @@ Node::Node(int ind) : Vector() @@
 Node::Node(const Vector &src) : Vector(src)
+{
   node_set = 0;
   index=(nnodes++);
   fixed=false;
   boundary=false;
   sam=false;
   dead = false;
+}
 Node::Node(double x,double y, double z) : Vector (x,y,z)
+{
   node_set = 0;
   index=(nnodes++);
   fixed=false;
   boundary=false;
   sam=false;
   dead = false;
+}
 Node::Node(const Node &src) : Vector(src)
+{
   node_set=0;
   owners=src.owners;
   m=src.m;
   index=src.index;
   fixed = src.fixed;
   boundary = src.boundary;
   sam=src.sam;
   dead = src.dead;
+}
 Cell &Node::getCell(const Neighbor &i)
+{
   return *i.getCell(); // use accessor!
+}
 ostream &Node::print(ostream &os) const {
   if (!dead)
     os << "Node ";
   else
     os << "DEAD NODE ";
   os  << index << "[ {" << x << ", " << y << ", " << z << "}: {";
   os << "Neighbors = { ";
   for (list<Neighbor>::const_iterator i =  owners.begin(); i!=owners.end(); i++) {
     os << " {" << i->cell->Index() << " " << i->nb1->Index() << " " << i->nb2->Index() << "} ";
+  }
   os << " } " << endl;
   return os;
+}
 #ifdef QTGRAPHICS
 void Node::DrawIndex(QGraphicsScene *c) const {
   //return DrawOwners(c);
   stringstream text;
   text << index;
   QGraphicsSimpleTextItem *number = new QGraphicsSimpleTextItem ( QString (text.str().c_str()), 0, c );
   number->setFont( QFont( "Helvetica", par.nodenumsize, QFont::Bold) );
   number->setPen( QPen (par.textcolor) );
   number->setZValue(20);
   number->show();
   Vector offs=Cell::Offset();
   number -> setPos(
 		   ((offs.x+x)*Cell::Factor()),
 		   ((offs.y+y)*Cell::Factor()) );
+}
 void Node::DrawOwners(QGraphicsScene *c) const {
   stringstream text;
   text << owners.size();
   QGraphicsSimpleTextItem *number = new QGraphicsSimpleTextItem ( QString (text.str().c_str()), 0, c );
   number->setFont( QFont( "Helvetica", par.nodenumsize, QFont::Bold) );
   number->setPen( QPen(par.textcolor) );
   number->setZValue(20);
   number->show();
   Vector offs=Cell::Offset();
   number ->setPos(((offs.x+x)*Cell::Factor()),
 		  ((offs.y+y)*Cell::Factor()) );
+}
 QVector<qreal> Node::NeighbourAngles(void)
+{
   QVector<qreal> angles;
   for (list<Neighbor>::iterator i=owners.begin();
        i!=owners.end();
        i++) {
   for (list<Neighbor>::iterator i=owners.begin(); i!=owners.end(); i++) {
     Vector v1 = (*this - *i->nb1).Normalised();
     Vector v2 = (*this - *i->nb2).Normalised();
     double angle = v1.SignedAngle(v2);
     if (angle<0) {
       //cerr << "Changing sign of angle " << angle << endl;
       angle = angle + 2*Pi;
+    }
     angles.push_back(angle);
     //cerr << "Cell " << i->cell->Index() << ": " <<  v1 << " and " << v2
     //     << ": angle = " << angles.back() << ", " << v1.Angle(v2) << endl;
+  }
   double sum=0.;
   for (QVector<qreal>::iterator i=angles.begin();
        i!=angles.end();
        i++) {
   for (QVector<qreal>::iterator i=angles.begin(); i!=angles.end(); i++) {
     sum+=*i;
+  }
   //cerr << "Angle sum = " << sum << endl;
   // Check if the summed angle is different from 2 Pi
   if (fabs(sum-(2*Pi))>0.0001) {
     MyWarning::warning("sum = %f",sum);
+  }
   return angles;
+}
 #endif
 ostream &operator<<(ostream &os, const Node &n) {
   n.print(os);
   return os;
+}
 /* finis */

src/parameter.cpp

➞

Show inline comments

@@ @@ -1787,96 +1787,97 @@ void Parameter::AssignValToPar(const cha @@
     if (s1) { free(s1); }
     s1=strdup(valc);
+  }
   if (!strcmp(namec, "s2")) {
     if (s2) { free(s2); }
     s2=strdup(valc);
+  }
   if (!strcmp(namec, "s3")) {
     if (s3) { free(s3); }
     s3=strdup(valc);
+  }
   if (!strcmp(namec, "b1")) {
     b1 = strtobool(valc);
+  }
   if (!strcmp(namec, "b2")) {
     b2 = strtobool(valc);
+  }
   if (!strcmp(namec, "b3")) {
     b3 = strtobool(valc);
+  }
   if (!strcmp(namec, "b4")) {
     b4 = strtobool(valc);
+  }
   if (!strcmp(namec, "dir1")) {
     if (dir1) { free(dir1); }
     dir1=strdup(valc);
+  }
   if (!strcmp(namec, "dir2")) {
     if (dir2) { free(dir2); }
     dir2=strdup(valc);
+  }
+}
 void Parameter::AssignValArrayToPar(const char *namec, vector<double> valarray) {
   if (!strcmp(namec, "D")) {
     int i=0;
     vector<double>::const_iterator v=valarray.begin();
     while (v!=valarray.end() && i <= 14 ) {
       D[i++]=*(v++);
+    }
+  }
   if (!strcmp(namec, "initval")) {
     int i=0;
     vector<double>::const_iterator v=valarray.begin();
     while (v!=valarray.end() && i <= 14 ) {
       initval[i++]=*(v++);
+    }
+  }
   if (!strcmp(namec, "k")) {
     int i=0;
     vector<double>::const_iterator v=valarray.begin();
     while (v!=valarray.end() && i <= 14 ) {
       k[i++]=*(v++);
+    }
+  }
+}
 /* void Parameter::XMLRead(xmlNode *root) {
    xmlNode *cur = root->xmlChildrenNode;
    while (cur!=NULL) {
    if ((!xmlStrcmp(cur->name, (const xmlChar *)"parameter"))){
    xmlNode *par_node = cur->xmlChildrenNode;
    while (par_node!=NULL) {
+   {
    if (!xmlStrcmp(par_node->name, (const xmlChar *)"par")) {
    xmlChar *namec = xmlGetProp(par_node, BAD_CAST "name");
    xmlChar *valc = xmlGetProp(par_node, BAD_CAST "val");
    if (valc) {
    AssignValToPar((const char*)namec,(const char*)valc);
    } else {
    // Probably a valarray
    xmlNode *sub_par_node = par_node->xmlChildrenNode;
    vector<double> valarray;
    while (sub_par_node != NULL) {
    if (!xmlStrcmp(sub_par_node->name, (const xmlChar *)"valarray")) {
    valarray = XMLIO::XMLReadValArray(sub_par_node);
+   }
    sub_par_node = sub_par_node->next;
+   }
    AssignValArrayToPar((const char*)namec, valarray);
+   }
+   }
+   }
    par_node = par_node->next;
+   }
+   }
    cur=cur->next;
+   }
    }*/
 ostream &operator<<(ostream &os, Parameter &p) {
   p.Write(os);
   return os;
+}
 /* finis */

src/parameter.h

➞

Show inline comments

@@ @@ -52,98 +52,99 @@ class Parameter { @@
   int resize_stride;
   double T;
   double lambda_length;
   double lambda_celllength;
   double target_length;
   double cell_expansion_rate;
   double cell_div_expansion_rate;
   bool auxin_dependent_growth;
   double ode_accuracy;
   double mc_stepsize;
   double mc_cell_stepsize;
   double energy_threshold;
   double bend_lambda;
   double alignment_lambda;
   double rel_cell_div_threshold;
   double rel_perimeter_stiffness;
   double collapse_node_threshold;
   double morphogen_div_threshold;
   double morphogen_expansion_threshold;
   bool copy_wall;
   double source;
   double * D;
   double * initval;
   double k1;
   double k2;
   double r;
   double kr;
   double km;
   double Pi_tot;
   double transport;
   double ka;
   double pin_prod;
   double pin_prod_in_epidermis;
   double pin_breakdown;
   double pin_breakdown_internal;
   double aux1prod;
   double aux1prodmeso;
   double aux1decay;
   double aux1decaymeso;
   double aux1transport;
   double aux_cons;
   double aux_breakdown;
   double kaux1;
   double kap;
   double leaf_tip_source;
   double sam_efflux;
   double sam_auxin;
   double sam_auxin_breakdown;
   double van3prod;
   double van3autokat;
   double van3sat;
   double k2van3;
   double dt;
   double rd_dt;
   char * datadir;
   bool movie;
   int nit;
   double maxt;
   int storage_stride;
   int xml_storage_stride;
   int rseed;
   int constituous_expansion_limit;
   double vessel_inh_level;
   double vessel_expansion_rate;
   double d;
   double e;
   double f;
   double c;
   double mu;
   double nu;
   double rho0;
   double rho1;
   double c0;
   double gamma;
   double eps;
   double * k;
   int i1;
   int i2;
   int i3;
   int i4;
   int i5;
   char * s1;
   char * s2;
   char * s3;
   bool b1;
   bool b2;
   bool b3;
   bool b4;
   char * dir1;
   char * dir2;
  private:
 };
 ostream &operator<<(ostream &os, Parameter &p);
 const char *sbool(const bool &p);
 #endif
-#endif
+/* finis */

src/pardialog.cpp

➞

Show inline comments

@@ @@ -710,96 +710,98 @@ void ParameterDialog::Reset(void) { @@
   lambda_length_edit->setText( QString("%1").arg(par.lambda_length) );
   lambda_celllength_edit->setText( QString("%1").arg(par.lambda_celllength) );
   target_length_edit->setText( QString("%1").arg(par.target_length) );
   cell_expansion_rate_edit->setText( QString("%1").arg(par.cell_expansion_rate) );
   cell_div_expansion_rate_edit->setText( QString("%1").arg(par.cell_div_expansion_rate) );
   auxin_dependent_growth_edit->setText( QString("%1").arg(sbool(par.auxin_dependent_growth)));
   ode_accuracy_edit->setText( QString("%1").arg(par.ode_accuracy) );
   mc_stepsize_edit->setText( QString("%1").arg(par.mc_stepsize) );
   mc_cell_stepsize_edit->setText( QString("%1").arg(par.mc_cell_stepsize) );
   energy_threshold_edit->setText( QString("%1").arg(par.energy_threshold) );
   bend_lambda_edit->setText( QString("%1").arg(par.bend_lambda) );
   alignment_lambda_edit->setText( QString("%1").arg(par.alignment_lambda) );
   rel_cell_div_threshold_edit->setText( QString("%1").arg(par.rel_cell_div_threshold) );
   rel_perimeter_stiffness_edit->setText( QString("%1").arg(par.rel_perimeter_stiffness) );
   collapse_node_threshold_edit->setText( QString("%1").arg(par.collapse_node_threshold) );
   morphogen_div_threshold_edit->setText( QString("%1").arg(par.morphogen_div_threshold) );
   morphogen_expansion_threshold_edit->setText( QString("%1").arg(par.morphogen_expansion_threshold) );
   copy_wall_edit->setText( QString("%1").arg(sbool(par.copy_wall)));
   source_edit->setText( QString("%1").arg(par.source) );
   QString D_string("%1,%2,%3,%4,%5,%6,%7,%8,%9,%10,%11,%12,%13,%14,%15");
   D_string = D_string.arg(par.D[0]).arg(par.D[1]).arg(par.D[2]).arg(par.D[3]).arg(par.D[4]).arg(par.D[5]).arg(par.D[6]).arg(par.D[7]).arg(par.D[8]).arg(par.D[9]).arg(par.D[10]).arg(par.D[11]).arg(par.D[12]).arg(par.D[13]).arg(par.D[14]);
   D_edit->setText( D_string );
   QString initval_string("%1,%2,%3,%4,%5,%6,%7,%8,%9,%10,%11,%12,%13,%14,%15");
   initval_string = initval_string.arg(par.initval[0]).arg(par.initval[1]).arg(par.initval[2]).arg(par.initval[3]).arg(par.initval[4]).arg(par.initval[5]).arg(par.initval[6]).arg(par.initval[7]).arg(par.initval[8]).arg(par.initval[9]).arg(par.initval[10]).arg(par.initval[11]).arg(par.initval[12]).arg(par.initval[13]).arg(par.initval[14]);
   initval_edit->setText( initval_string );
   k1_edit->setText( QString("%1").arg(par.k1) );
   k2_edit->setText( QString("%1").arg(par.k2) );
   r_edit->setText( QString("%1").arg(par.r) );
   kr_edit->setText( QString("%1").arg(par.kr) );
   km_edit->setText( QString("%1").arg(par.km) );
   Pi_tot_edit->setText( QString("%1").arg(par.Pi_tot) );
   transport_edit->setText( QString("%1").arg(par.transport) );
   ka_edit->setText( QString("%1").arg(par.ka) );
   pin_prod_edit->setText( QString("%1").arg(par.pin_prod) );
   pin_prod_in_epidermis_edit->setText( QString("%1").arg(par.pin_prod_in_epidermis) );
   pin_breakdown_edit->setText( QString("%1").arg(par.pin_breakdown) );
   pin_breakdown_internal_edit->setText( QString("%1").arg(par.pin_breakdown_internal) );
   aux1prod_edit->setText( QString("%1").arg(par.aux1prod) );
   aux1prodmeso_edit->setText( QString("%1").arg(par.aux1prodmeso) );
   aux1decay_edit->setText( QString("%1").arg(par.aux1decay) );
   aux1decaymeso_edit->setText( QString("%1").arg(par.aux1decaymeso) );
   aux1transport_edit->setText( QString("%1").arg(par.aux1transport) );
   aux_cons_edit->setText( QString("%1").arg(par.aux_cons) );
   aux_breakdown_edit->setText( QString("%1").arg(par.aux_breakdown) );
   kaux1_edit->setText( QString("%1").arg(par.kaux1) );
   kap_edit->setText( QString("%1").arg(par.kap) );
   leaf_tip_source_edit->setText( QString("%1").arg(par.leaf_tip_source) );
   sam_efflux_edit->setText( QString("%1").arg(par.sam_efflux) );
   sam_auxin_edit->setText( QString("%1").arg(par.sam_auxin) );
   sam_auxin_breakdown_edit->setText( QString("%1").arg(par.sam_auxin_breakdown) );
   van3prod_edit->setText( QString("%1").arg(par.van3prod) );
   van3autokat_edit->setText( QString("%1").arg(par.van3autokat) );
   van3sat_edit->setText( QString("%1").arg(par.van3sat) );
   k2van3_edit->setText( QString("%1").arg(par.k2van3) );
   dt_edit->setText( QString("%1").arg(par.dt) );
   rd_dt_edit->setText( QString("%1").arg(par.rd_dt) );
   datadir_edit->setText( QString("%1").arg(par.datadir) );
   movie_edit->setText( QString("%1").arg(sbool(par.movie)));
   nit_edit->setText( QString("%1").arg(par.nit) );
   maxt_edit->setText( QString("%1").arg(par.maxt) );
   storage_stride_edit->setText( QString("%1").arg(par.storage_stride) );
   xml_storage_stride_edit->setText( QString("%1").arg(par.xml_storage_stride) );
   rseed_edit->setText( QString("%1").arg(par.rseed) );
   constituous_expansion_limit_edit->setText( QString("%1").arg(par.constituous_expansion_limit) );
   vessel_inh_level_edit->setText( QString("%1").arg(par.vessel_inh_level) );
   vessel_expansion_rate_edit->setText( QString("%1").arg(par.vessel_expansion_rate) );
   d_edit->setText( QString("%1").arg(par.d) );
   e_edit->setText( QString("%1").arg(par.e) );
   f_edit->setText( QString("%1").arg(par.f) );
   c_edit->setText( QString("%1").arg(par.c) );
   mu_edit->setText( QString("%1").arg(par.mu) );
   nu_edit->setText( QString("%1").arg(par.nu) );
   rho0_edit->setText( QString("%1").arg(par.rho0) );
   rho1_edit->setText( QString("%1").arg(par.rho1) );
   c0_edit->setText( QString("%1").arg(par.c0) );
   gamma_edit->setText( QString("%1").arg(par.gamma) );
   eps_edit->setText( QString("%1").arg(par.eps) );
   QString k_string("%1,%2,%3,%4,%5,%6,%7,%8,%9,%10,%11,%12,%13,%14,%15");
   k_string = k_string.arg(par.k[0]).arg(par.k[1]).arg(par.k[2]).arg(par.k[3]).arg(par.k[4]).arg(par.k[5]).arg(par.k[6]).arg(par.k[7]).arg(par.k[8]).arg(par.k[9]).arg(par.k[10]).arg(par.k[11]).arg(par.k[12]).arg(par.k[13]).arg(par.k[14]);
   k_edit->setText( k_string );
   i1_edit->setText( QString("%1").arg(par.i1) );
   i2_edit->setText( QString("%1").arg(par.i2) );
   i3_edit->setText( QString("%1").arg(par.i3) );
   i4_edit->setText( QString("%1").arg(par.i4) );
   i5_edit->setText( QString("%1").arg(par.i5) );
   s1_edit->setText( QString("%1").arg(par.s1) );
   s2_edit->setText( QString("%1").arg(par.s2) );
   s3_edit->setText( QString("%1").arg(par.s3) );
   b1_edit->setText( QString("%1").arg(sbool(par.b1)));
   b2_edit->setText( QString("%1").arg(sbool(par.b2)));
   b3_edit->setText( QString("%1").arg(sbool(par.b3)));
   b4_edit->setText( QString("%1").arg(sbool(par.b4)));
   dir1_edit->setText( QString("%1").arg(par.dir1) );
   dir2_edit->setText( QString("%1").arg(par.dir2) );
+}
 /* finis */

src/pardialog.h

➞

Show inline comments

@@ @@ -49,97 +49,100 @@ class ParameterDialog : public QDialog { @@
   QLineEdit *cellnumsize_edit;
   QLineEdit *nodenumsize_edit;
   QLineEdit *node_mag_edit;
   QLineEdit *outlinewidth_edit;
   QLineEdit *cell_outline_color_edit;
   QLineEdit *resize_stride_edit;
   QLineEdit *T_edit;
   QLineEdit *lambda_length_edit;
   QLineEdit *lambda_celllength_edit;
   QLineEdit *target_length_edit;
   QLineEdit *cell_expansion_rate_edit;
   QLineEdit *cell_div_expansion_rate_edit;
   QLineEdit *auxin_dependent_growth_edit;
   QLineEdit *ode_accuracy_edit;
   QLineEdit *mc_stepsize_edit;
   QLineEdit *mc_cell_stepsize_edit;
   QLineEdit *energy_threshold_edit;
   QLineEdit *bend_lambda_edit;
   QLineEdit *alignment_lambda_edit;
   QLineEdit *rel_cell_div_threshold_edit;
   QLineEdit *rel_perimeter_stiffness_edit;
   QLineEdit *collapse_node_threshold_edit;
   QLineEdit *morphogen_div_threshold_edit;
   QLineEdit *morphogen_expansion_threshold_edit;
   QLineEdit *copy_wall_edit;
   QLineEdit *source_edit;
   QLineEdit *D_edit;
   QLineEdit *initval_edit;
   QLineEdit *k1_edit;
   QLineEdit *k2_edit;
   QLineEdit *r_edit;
   QLineEdit *kr_edit;
   QLineEdit *km_edit;
   QLineEdit *Pi_tot_edit;
   QLineEdit *transport_edit;
   QLineEdit *ka_edit;
   QLineEdit *pin_prod_edit;
   QLineEdit *pin_prod_in_epidermis_edit;
   QLineEdit *pin_breakdown_edit;
   QLineEdit *pin_breakdown_internal_edit;
   QLineEdit *aux1prod_edit;
   QLineEdit *aux1prodmeso_edit;
   QLineEdit *aux1decay_edit;
   QLineEdit *aux1decaymeso_edit;
   QLineEdit *aux1transport_edit;
   QLineEdit *aux_cons_edit;
   QLineEdit *aux_breakdown_edit;
   QLineEdit *kaux1_edit;
   QLineEdit *kap_edit;
   QLineEdit *leaf_tip_source_edit;
   QLineEdit *sam_efflux_edit;
   QLineEdit *sam_auxin_edit;
   QLineEdit *sam_auxin_breakdown_edit;
   QLineEdit *van3prod_edit;
   QLineEdit *van3autokat_edit;
   QLineEdit *van3sat_edit;
   QLineEdit *k2van3_edit;
   QLineEdit *dt_edit;
   QLineEdit *rd_dt_edit;
   QLineEdit *datadir_edit;
   QLineEdit *movie_edit;
   QLineEdit *nit_edit;
   QLineEdit *maxt_edit;
   QLineEdit *storage_stride_edit;
   QLineEdit *xml_storage_stride_edit;
   QLineEdit *rseed_edit;
   QLineEdit *constituous_expansion_limit_edit;
   QLineEdit *vessel_inh_level_edit;
   QLineEdit *vessel_expansion_rate_edit;
   QLineEdit *d_edit;
   QLineEdit *e_edit;
   QLineEdit *f_edit;
   QLineEdit *c_edit;
   QLineEdit *mu_edit;
   QLineEdit *nu_edit;
   QLineEdit *rho0_edit;
   QLineEdit *rho1_edit;
   QLineEdit *c0_edit;
   QLineEdit *gamma_edit;
   QLineEdit *eps_edit;
   QLineEdit *k_edit;
   QLineEdit *i1_edit;
   QLineEdit *i2_edit;
   QLineEdit *i3_edit;
   QLineEdit *i4_edit;
   QLineEdit *i5_edit;
   QLineEdit *s1_edit;
   QLineEdit *s2_edit;
   QLineEdit *s3_edit;
   QLineEdit *b1_edit;
   QLineEdit *b2_edit;
   QLineEdit *b3_edit;
   QLineEdit *b4_edit;
   QLineEdit *dir1_edit;
   QLineEdit *dir2_edit;
 };
 #endif
 /* finis */

src/xmlwrite.cpp

➞

Show inline comments

@@ @@ -657,553 +657,539 @@ void Wall::XMLAdd(xmlNode *parent) const @@
   if (apoplast) {
     xmlNodePtr apo_xml = xmlNewChild(xmlwall, NULL, BAD_CAST "apoplast", NULL);
     for (int i=0;i<Cell::NChem();i++) {
       xmlNodePtr apo_val_xml = xmlNewChild(apo_xml, NULL, BAD_CAST "val", NULL);
+      {
 	ostringstream text;
 	text << transporters2[i];
 	xmlNewProp(apo_val_xml, BAD_CAST "v", BAD_CAST text.str().c_str());
+      }
+    }
+  }
+}
 vector<double> XMLIO::XMLReadValArray(xmlNode *cur)
+{
   vector<double> result;
   QLocale standardlocale(QLocale::C);
   bool ok;
   xmlNode *valarray_node = cur->xmlChildrenNode;
   while (valarray_node!=NULL) {
     if (!xmlStrcmp(valarray_node->name, (const xmlChar *)"val")) {
       xmlChar *vc = xmlGetProp(valarray_node, BAD_CAST "v");
       if (vc) {
 	//result.push_back(strtod( (const char *)vc, 0));
 	result.push_back(standardlocale.toDouble((char *)vc, &ok));
 	if (!ok) MyWarning::error("Could not convert \"%s\" to double in XMLRead.",(char *)vc);
 	xmlFree(vc);
+      }
+    }
     valarray_node = valarray_node->next;
+  }
   return result;
+}
 void Mesh::XMLSave(const char *docname, xmlNode *options) const
+{
   // based on libxml2 example code "tree2.c"
   xmlDocPtr doc = NULL;       /* document pointer */
   xmlNodePtr root_node = NULL;/* node pointers */
   //xmlDtdPtr dtd = NULL;       /* DTD pointer */
   //  LIBXML_TEST_VERSION;
   /*
    * Creates a new document, a node and set it as a root node
    */
   doc = xmlNewDoc(BAD_CAST "1.0");
   root_node = xmlNewNode(NULL, BAD_CAST "leaf");
   xmlDocSetRootElement(doc, root_node);
   /*
    * xmlNewProp() creates attributes, which is "attached" to an node.
    * It returns xmlAttrPtr, which isn't used here.
    */
   xmlNewProp(root_node, BAD_CAST "name", BAD_CAST docname);
   time_t t;
   std::time(&t);
   char *tstring = strdup(asctime(localtime(&t))); // but this does
   // replace "end of line character by '\0'
   char *eol=strchr(tstring,'\n');
   if (eol!=NULL)
     *eol='\0';
   xmlNewProp(root_node, BAD_CAST "date", BAD_CAST tstring);
   free(tstring);
   QString simtime = QString("%1").arg(time);
   xmlNewProp(root_node, BAD_CAST "simtime", BAD_CAST simtime.toStdString().c_str());
   /*
    * Creates a DTD declaration. Isn't mandatory.
    */
   par.XMLAdd(root_node);
   xmlNodePtr xmlnodes = xmlNewChild(root_node, NULL, BAD_CAST "nodes",NULL);
   { ostringstream text;
     text << NNodes();
     xmlNewProp(xmlnodes, BAD_CAST "n", BAD_CAST text.str().c_str());
+  }
   { ostringstream text;
     text << Node::target_length;
     xmlNewProp(xmlnodes, BAD_CAST "target_length", BAD_CAST text.str().c_str());
+  }
   for (vector<Node *>::const_iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (vector<Node *>::const_iterator i=nodes.begin(); i!=nodes.end(); i++) {
     (*i)->XMLAdd(xmlnodes) ;
+  }
   /*
    * xmlNewChild() creates a new node, which is "attached" as child node
    * of root_node node.
    */
   xmlNodePtr xmlcells = xmlNewChild(root_node, NULL, BAD_CAST "cells",NULL);
+  {
     ostringstream text;
     text << NCells();
     xmlNewProp(xmlcells, BAD_CAST "n", BAD_CAST text.str().c_str());
+  }
+  {
     ostringstream text;
     text << Cell::offset[0];
     xmlNewProp(xmlcells, BAD_CAST "offsetx", BAD_CAST text.str().c_str());
+  }
+  {
     ostringstream text;
     text << Cell::offset[1];
     xmlNewProp(xmlcells, BAD_CAST "offsety", BAD_CAST text.str().c_str());
+  }
+  {
     ostringstream text;
     text << Cell::factor;
     xmlNewProp(xmlcells, BAD_CAST "magnification", BAD_CAST text.str().c_str());
+  }
+  {
     ostringstream text;
     text << cells.front()->BaseArea();
     xmlNewProp(xmlcells, BAD_CAST "base_area", BAD_CAST text.str().c_str());
+  }
+  {
     ostringstream text;
     text << Cell::NChem();
     xmlNewProp(xmlcells, BAD_CAST "nchem", BAD_CAST text.str().c_str());
+  }
   for (vector<Cell *>::const_iterator i=cells.begin();
        i!=cells.end();
        i++) {
   for (vector<Cell *>::const_iterator i=cells.begin(); i!=cells.end(); i++) {
     (*i)->XMLAdd(xmlcells) ;
+  }
   boundary_polygon->XMLAdd(xmlcells);
   xmlNodePtr xmlwalls = xmlNewChild(root_node, NULL, BAD_CAST "walls",NULL);
+  {
     ostringstream text;
     text << walls.size();
     xmlNewProp(xmlwalls, BAD_CAST "n", BAD_CAST text.str().c_str());
+  }
   for (list<Wall *>::const_iterator i=walls.begin();
        i!=walls.end();
        i++) {
   for (list<Wall *>::const_iterator i=walls.begin(); i!=walls.end(); i++) {
     (*i)->XMLAdd(xmlwalls) ;
+  }
   xmlNodePtr xmlnodesets = xmlNewChild(root_node, NULL, BAD_CAST "nodesets",NULL);
+  {
     ostringstream text;
     text << node_sets.size();
     xmlNewProp(xmlnodesets, BAD_CAST "n", BAD_CAST text.str().c_str());
+  }
   for_each( node_sets.begin(), node_sets.end(), bind2nd ( mem_fun( &NodeSet::XMLAdd ), xmlnodesets ) );
   // Add option tree for interactive application
   if (options) {
     xmlAddChild(root_node, options);
+  }
   /*
    * Dumping document to stdio or file
    */
   xmlSetDocCompressMode(doc,9);
   xmlSaveFormatFileEnc(docname, doc, "UTF-8", 1);
   /*free the document */
   xmlFreeDoc(doc);
   /*
    *Free the global variables that may
    *have been allocated by the parser.
    */
   xmlCleanupParser();
   /*
    * this is to debug memory for regression tests
    */
+}
 void Mesh::XMLReadSimtime(const xmlNode *a_node)
+{
   xmlNode *root_node;
   root_node = (xmlNode *)a_node;
   xmlChar *strsimtime = xmlGetProp(root_node, BAD_CAST "simtime");
   if (strsimtime) {
     QLocale standardlocale(QLocale::C);
     bool ok;
     double simtime=standardlocale.toDouble((char *)strsimtime, &ok);
     if (!ok) MyWarning::error("Could not convert \"%s\" to double in XMLRead.",(char *)strsimtime);
     time = simtime;
 #ifdef QDEBUG
     qDebug() << "Simtime = " << strsimtime << endl;
 #endif
   } else {
 #ifdef QDEBUG
     qDebug() << "No simtime found in file." << endl;
 #endif
     time =0;
+  }
+}
 void Mesh::XMLReadPars(const xmlNode * a_node)
+{
   xmlNode *root_node;
   root_node = (xmlNode *)a_node;
   par.XMLRead(root_node);
   Seed(par.rseed);
   MakeDir(par.datadir);
+}
 void Mesh::XMLReadGeometry(const xmlNode * a_node)
+{
   xmlNode *root_node;
   root_node = (xmlNode *)a_node;
   xmlNode *cur;
   QLocale standardlocale(QLocale::C);
   bool ok;
   // allocate Nodes
   cur = root_node->xmlChildrenNode;
   while (cur!=NULL) {
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"nodes"))){
       XMLReadNodes(cur);
+    }
     cur=cur->next;
+  }
   cur = root_node->xmlChildrenNode;
   while (cur!=NULL) {
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"nodesets"))) {
       XMLReadNodeSets(cur);
+    }
     cur=cur->next;
+  }
   cur=root_node->xmlChildrenNode;
   while (cur!=NULL) {
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"nodes"))) {
       XMLReadNodeSetsToNodes(cur);
+    }
     cur = cur->next;
+  }
   // allocate Cells
   cur = root_node;
   // allocate Cells
   cur = cur->xmlChildrenNode;
   while (cur!=NULL) {
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"cells"))){
       xmlChar *offsetxc = xmlGetProp(cur, BAD_CAST "offsetx");
       xmlChar *offsetyc = xmlGetProp(cur, BAD_CAST "offsety");
       double ox=standardlocale.toDouble((const char *)offsetxc, &ok);
       if (!ok) MyWarning::error("Could not convert \"%s\" to double in XMLRead.",(const char *)offsetxc);
       double oy=standardlocale.toDouble((const char *)offsetyc, &ok);
       if (!ok) MyWarning::error("Could not convert \"%s\" to double in XMLRead.",(const char *)offsetyc);
       Cell::setOffset(ox, oy);
       xmlFree(offsetxc);
       xmlFree(offsetyc);
       xmlChar *magnificationc = xmlGetProp(cur, BAD_CAST "magnification");
       Cell::SetMagnification(standardlocale.toDouble((const char *)magnificationc, &ok));
       if (!ok) MyWarning::error("Could not convert \"%s\" to double in XMLRead.",(const char *)magnificationc);
       xmlFree(magnificationc);
       xmlChar *baseareac = xmlGetProp(cur, BAD_CAST "base_area");
       Cell::BaseArea() = standardlocale.toDouble((const char *)baseareac, &ok);
       if (!ok) MyWarning::error("Could not convert \"%s\" to double in XMLRead.",(const char *)baseareac);
       xmlFree(baseareac);
       XMLReadCells(cur);
+    }
     cur=cur->next;
+  }
   // allocate Walls (we need to have read the cells before constructing walls)
   vector <Wall *> tmp_walls;
   cur = root_node->xmlChildrenNode;
   while (cur!=NULL) {
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"walls"))){
       XMLReadWalls(cur, &tmp_walls);
+    }
     cur=cur->next;
+  }
   // read walls to cells and boundary_polygon
   cur = root_node->xmlChildrenNode;
   while (cur!=NULL) {
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"cells"))) {
       XMLReadWallsToCells(cur, &tmp_walls);
+    }
     cur=cur->next;
+  }
   boundary_polygon->ConstructNeighborList();
   boundary_polygon->ConstructConnections();
   for (vector<Cell *>::iterator c=cells.begin();
        c!=cells.end();
        c++) {
   for (vector<Cell *>::iterator c=cells.begin(); c!=cells.end(); c++) {
     (*c)->ConstructNeighborList();
     (*c)->ConstructConnections();
+  }
   shuffled_cells.clear();
   shuffled_cells = cells;
   MyUrand r(shuffled_cells.size());
   random_shuffle(shuffled_cells.begin(),shuffled_cells.end(),r);
+}
 void Mesh::XMLParseTree(const xmlNode *root_node)
+{
   XMLReadSimtime(root_node);
   XMLReadPars(root_node);
   XMLReadGeometry(root_node);
+}
 void Mesh::XMLReadNodes(xmlNode *root)
+{
   QLocale standardlocale(QLocale::C);
   bool ok;
   xmlNode *cur = root;
   cur = cur->xmlChildrenNode;
   for (vector<Node *>::iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
   for (vector<Node *>::iterator i=nodes.begin(); i!=nodes.end(); i++) {
     delete *i;
+  }
   nodes.clear();
   Node::nnodes=0;
   xmlChar *tlc = xmlGetProp(root, BAD_CAST "target_length");
   if (tlc != 0) {
     Node::target_length = standardlocale.toDouble((const char *)tlc, &ok);
     if (!ok) MyWarning::error("Could not convert \"%s\" to double in XMLRead.",(const char *)tlc);
     xmlFree(tlc);
   } else {
     // note that libxml2 also defines a token "warning"
     MyWarning::unique_warning("Warning: value found in XML file for Node::target_length.");
+  }
   while (cur!=NULL) {
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"node"))){
       xmlChar *xc = xmlGetProp(cur, BAD_CAST "x");
       if (xc==0) {
 	unique_warning("Token \"x\" not found in xmlwrite.cpp at or around line no. 722");
+      }
       xmlChar *yc = xmlGetProp(cur, BAD_CAST "y");
       if (yc==0) {
 	unique_warning("Token \"y\" not found in xmlwrite.cpp at or around line no. 727");
+      }
       xmlChar *fixedc = xmlGetProp(cur, BAD_CAST "fixed");
       if (fixedc==0) {
 	unique_warning("Token \"fixed\" not found in xmlwrite.cpp at or around line.");
+      }
       xmlChar *boundaryc = xmlGetProp(cur, BAD_CAST "boundary");
       if (boundaryc==0) {
 	unique_warning("Token \"boundary\" not found in xmlwrite.cpp at or around line.");
+      }
       xmlChar *samc = xmlGetProp(cur, BAD_CAST "sam");
       if (samc==0) {
 	unique_warning("Token \"sam\" not found in xmlwrite.cpp at or around line.");
+      }
       double x = standardlocale.toDouble((const char *)xc, &ok);
       if (!ok) MyWarning::error("Could not convert \"%s\" to double in XMLRead.",(const char *)xc);
       double y = standardlocale.toDouble((const char *)yc, &ok);
       if (!ok) MyWarning::error("Could not convert \"%s\" to double in XMLRead.",(const char *)yc);
       Node *new_node = new Node(x,y);
       nodes.push_back(new_node);
       new_node->m = this;
       new_node->fixed = strtobool( (char *)fixedc);
       new_node->boundary = strtobool( (char *)boundaryc );
       new_node->sam = strtobool ( (char *)samc);
       new_node->node_set = 0;
       xmlFree(xc);
       xmlFree(yc);
       xmlFree(boundaryc);
       xmlFree(fixedc);
       xmlFree(samc);
+    }
     cur=cur->next;
+  }
   shuffled_nodes.clear();
   shuffled_nodes = nodes;
   MyUrand r(shuffled_nodes.size());
   random_shuffle(shuffled_nodes.begin(),shuffled_nodes.end(),r);
+}
 void Mesh::XMLReadWalls(xmlNode *root, vector<Wall *> *tmp_walls)
+{
   xmlNode *cur = root;
   cur = cur->xmlChildrenNode;
   for (list<Wall *>::iterator i=walls.begin();
        i!=walls.end();
        i++) {
   for (list<Wall *>::iterator i=walls.begin(); i!=walls.end(); i++) {
     delete *i;
+  }
   walls.clear();
   Wall::nwalls = 0;
   tmp_walls->clear();
   QLocale standardlocale(QLocale::C);
   bool ok;
   while (cur!=NULL) {
     vector<int> tmp_nodes;
     while(cur!=NULL) {
       if ((!xmlStrcmp(cur->name, (const xmlChar *)"wall"))) {
 	xmlChar *nc = xmlGetProp(cur, BAD_CAST "c1");
 	if (nc==0) {
 	  unique_warning("Token \"c1\" not found in xmlwrite.cpp at or around line no. 777");
+	}
 	int c1 = atoi( (char *)nc);
 	xmlFree(nc);
 	nc = xmlGetProp(cur, BAD_CAST "c2");
 	if (	nc==0) {
 	  unique_warning("Token \"c2\" not found in xmlwrite.cpp at or around line no. 785");
+	}
 	int c2 = atoi( (char *)nc);
 	xmlFree(nc);
 	nc = xmlGetProp(cur, BAD_CAST "n1");
 	if (	nc==0) {
 	  unique_warning("Token \"n1\" not found in xmlwrite.cpp at or around line no. 793");
+	}
 	int n1 = atoi( (char *)nc);
 	xmlFree(nc);
 	nc = xmlGetProp(cur, BAD_CAST "n2");
 	if (	nc==0) {
 	  unique_warning("Token \"n2\" not found in xmlwrite.cpp at or around line no. 801");
+	}
 	int n2 = atoi( (char *)nc);
 	xmlFree(nc);
 	nc = xmlGetProp(cur, BAD_CAST "length");
 	if (	nc==0) {
 	  unique_warning("Token \"length\" not found in xmlwrite.cpp at or around line no. 809");
+	}
 	double length = standardlocale.toDouble((const char *)nc, &ok);
 	if (!ok) MyWarning::error("Could not convert \"%s\" to double in XMLRead.",(const char *)nc);
 	xmlFree(nc);
 	nc = xmlGetProp(cur, BAD_CAST "viz_flux");
 	double viz_flux = 0.0;
 	if (nc!=0) {
 	  viz_flux = standardlocale.toDouble((const char *)nc, &ok);
 	  if (!ok) MyWarning::error("Could not convert \"%s\" to double in XMLRead.",(const char *)nc);
+	}
 	xmlFree(nc);
 	Wall::WallType wall_type;
+	{
 	  xmlChar *v_str = xmlGetProp(cur, BAD_CAST "wall_type");
 	  if (v_str != 0) {
 	    if (!xmlStrcmp(v_str, (const xmlChar *)"aux_source")) {
 	      wall_type = Wall::AuxSource;
 	    } else {
 	      if (!xmlStrcmp(v_str, (const xmlChar *)"aux_sink")) {
 		wall_type = Wall::AuxSink;
 	      } else {
 		wall_type = Wall::Normal;
+	      }
+	    }
 	    xmlFree(v_str);
 	  } else {
 	    wall_type = Wall::Normal;
+	  }
+	}
 	bool dead = false;
+	{
 	  // Note: property "delete" is used to manually clean up wall lists in XML files
 	  // Simply add property 'delete="true"' to the wall and it will be removed from
@@ @@ -1332,235 +1318,229 @@ void Mesh::XMLReadWalls(xmlNode *root, v @@
+  }
+}
 void Mesh::XMLReadWallsToCells(xmlNode *root, vector<Wall *> *tmp_walls)
+{
   // Add the walls to the cells (do this after reading the walls; read walls after reading cells...)
   // 1. Read Nodes
   // 2. Read Cells
   // 3. Read Walls
   // 4. Read Walls into Cells
   xmlNode *cur = root->xmlChildrenNode;
   int ci=0; // cell index
   while (cur!=NULL) {
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"cell")) ||
 	(!xmlStrcmp(cur->name, (const xmlChar *)"boundary_polygon" ))) {
       vector<int> tmp_walls_ind;
       xmlNode *n = cur->xmlChildrenNode;
       while(n!=NULL) {
 	if ((!xmlStrcmp(n->name, (const xmlChar *)"wall"))) {
 	  xmlChar *nc = xmlGetProp(n, BAD_CAST "w");
 	  if (nc==0) {
 	    unique_warning("Token \"w\" not found in xmlwrite.cpp at or around line no. 931");
+	  }
 	  tmp_walls_ind.push_back(atoi( (char *)nc));
 	  xmlFree(nc);
+	}
 	n = n->next;
+      }
       if (!xmlStrcmp(cur->name, (const xmlChar *)"boundary_polygon")) {
 	int nwalls = tmp_walls_ind.size();
 	for (int i=0;i<nwalls;i++) {
 	  boundary_polygon->walls.push_back((*tmp_walls)[tmp_walls_ind[i]]);
+	}
       } else {
 	int nwalls = tmp_walls_ind.size();
 	for (int i=0;i<nwalls;i++) {
 	  cells[ci]->walls.push_back((*tmp_walls)[tmp_walls_ind[i]]);
+	}
 	ci++;
+      }
+    }
     cur=cur->next;
+  }
+}
 void Mesh::XMLReadNodeSetsToNodes(xmlNode *root)
+{
   xmlNode *cur = root->xmlChildrenNode;
   int ci=0; // cell index
   while (cur!=NULL) {
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"node"))) {
       xmlNode *n = cur->xmlChildrenNode;
       while(n!=NULL) {
 	xmlChar *nc = xmlGetProp(n, BAD_CAST "nodeset");
 	if (nc!=0) {
 	  int nodeset_n = atoi( (char *)nc);
 	  nodes[ci]->node_set = node_sets[nodeset_n];
 	  xmlFree(nc);
 	} else {
 	  nodes[ci]->node_set = 0;
+	}
 	n = n->next;
+      }
       ci++;
+    }
     cur=cur->next;
+  }
+}
 void Mesh::XMLReadNodeSets(xmlNode *root)
+{
   for (vector<NodeSet *>::iterator i=node_sets.begin();
        i!=node_sets.end();
        i++) {
   for (vector<NodeSet *>::iterator i=node_sets.begin(); i!=node_sets.end(); i++) {
     delete *i;
+  }
   node_sets.clear();
   xmlNode *cur = root->xmlChildrenNode;
   while (cur!=NULL) {
     NodeSet *new_nodeset=0;
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"nodeset"))){
       new_nodeset = new NodeSet();
       node_sets.push_back(new_nodeset);
+    }
     if (new_nodeset == 0) {
       cur = cur->next;
       continue;
+    }
     new_nodeset->XMLRead(cur, this);
     cur=cur->next;
+  }
+}
 void Mesh::XMLReadCells(xmlNode *root)
+{
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
   for (vector<Cell *>::iterator i=cells.begin(); i!=cells.end(); i++) {
     delete *i;
+  }
   cells.clear();
   Cell::NCells() = 0;
   delete boundary_polygon;
   xmlNode *cur = root->xmlChildrenNode;
   while (cur!=NULL) {
     Cell *new_cell=0;
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"cell"))){
       new_cell = new Cell(0,0);
       new_cell->m = this;
       cells.push_back(new_cell);
     } else {
       if ((!xmlStrcmp(cur->name, (const xmlChar *)"boundary_polygon"))) {
 	new_cell = boundary_polygon = new BoundaryPolygon(0,0);
 	boundary_polygon->m = this;
+      }
+    }
     if (new_cell == 0) {
       cur = cur->next;
       continue;
+    }
     new_cell->XMLRead(cur);
     cur=cur->next;
+  }
+}
 void Mesh::XMLRead(const char *docname, xmlNode **settings, bool geometry, bool pars, bool simtime)
+{
   xmlDocPtr doc = xmlParseFile(docname);
   if (doc == NULL ) {
     throw("Document not parsed successfully.");
     return;
+  }
   xmlNodePtr cur = xmlDocGetRootElement(doc);
   if (cur == NULL) {
     throw("Document is empty");
     xmlFreeDoc(doc);
     return;
+  }
   if (xmlStrcmp(cur->name, (const xmlChar *) "leaf")) {
     throw("XML file of the wrong type, it is not a leaf.");
     xmlFreeDoc(doc);
     return;
+  }
   /*Get the root element node */
   xmlNode *root_element = xmlDocGetRootElement(doc);
   if (geometry) XMLReadGeometry(root_element);
   if (pars) XMLReadPars(root_element);
   if (simtime) XMLReadSimtime(root_element);
   // If pointer settings defined, return a copy of the settings tree
   if (settings) {
     xmlNode *cur = root_element->xmlChildrenNode;
     // if settings field is not found, *settings will be returned 0.
     *settings = 0;
     while (cur!=NULL) {
       if ((!xmlStrcmp(cur->name, (const xmlChar *)"settings"))){
 	*settings = xmlCopyNode(cur,1);
+      }
       cur=cur->next;
+    }
+  }
   xmlFreeDoc(doc);
   /*
    *Free the global variables that may
    *have been allocated by the parser.
    */
   xmlCleanupParser();
   // We're doing this so we can manually delete walls with by adding the 'delete="true"' property
   CleanUpCellNodeLists();
+}
 void Parameter::XMLRead(xmlNode *root)
+{
   xmlNode *cur = root->xmlChildrenNode;
   while (cur!=NULL) {
     if ((!xmlStrcmp(cur->name, (const xmlChar *)"parameter"))){
       xmlNode *par_node = cur->xmlChildrenNode;
       while (par_node!=NULL) {
+	{
 	  if (!xmlStrcmp(par_node->name, (const xmlChar *)"par")) {
 	    xmlChar *namec = xmlGetProp(par_node, BAD_CAST "name");
 	    xmlChar *valc = xmlGetProp(par_node, BAD_CAST "val");
 	    if (valc) {

0 comments (0 inline, 0 general)