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

@@ @@ -76,51 +76,49 @@ public: @@
 };
 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);
+	}

src/canvas.cpp

➞

Show inline comments

@@ @@ -274,100 +274,92 @@ void FigureEditor::mouseReleaseEvent(QMo @@
       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

src/cell.cpp

➞

Show inline comments

@@ @@ -61,88 +61,84 @@ bool Cell::Eq(Cell *c) const { return th @@
 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
@@ @@ -152,205 +148,197 @@ void Cell::Apoptose(void) @@
 #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
@@ @@ -396,51 +384,49 @@ void Cell::DivideWalls(ItList new_node_l @@
 #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);
@@ @@ -673,51 +659,49 @@ void Cell::DivideWalls(ItList new_node_l @@
     // 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.
@@ @@ -956,267 +940,234 @@ void Cell::DivideWalls(ItList new_node_l @@
   // 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;
@@ @@ -1248,52 +1199,49 @@ bool Cell::MoveSelfIntersectsP(Node *mov @@
   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;
@@ @@ -1334,267 +1282,246 @@ bool Cell::IntersectsWithLineP(const Vec @@
 	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 {
@@ @@ -1662,115 +1589,99 @@ void Cell::DrawFluxes(QGraphicsScene *c, @@
 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

@@ @@ -213,179 +213,169 @@ CellBase CellBase::operator=(const CellB @@
   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)*
@@ @@ -431,51 +421,49 @@ double CellBase::Length(Vector *long_axi @@
     //   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)*
@@ @@ -582,69 +570,60 @@ void CellBase::Dump(ostream &os) const @@
+  }
   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

@@ @@ -44,96 +44,90 @@ @@
 #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;
@@ @@ -223,79 +217,49 @@ Cell &Mesh::EllipticCell(double xc, doub @@
+  }
   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));
@@ @@ -310,174 +274,140 @@ Cell &Mesh::LeafPrimordium(int nnodes, d @@
 		  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;
@@ @@ -491,128 +421,118 @@ void Mesh::SetBaseArea(void) { @@
   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 )) {
@@ @@ -621,51 +541,49 @@ double Mesh::DisplaceNodes(void) { @@
     }*/
     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);
@@ @@ -939,67 +857,61 @@ double Mesh::DisplaceNodes(void) { @@
 	//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++) ) {
@@ @@ -1214,294 +1126,224 @@ void Mesh::CircumCircle(double x1,double @@
     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;
+  }
 };
@@ @@ -1575,51 +1417,49 @@ void Mesh::RepairBoundaryPolygon(void) { @@
     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;
@@ @@ -1681,167 +1521,144 @@ Node* Mesh::findNextBoundaryNode(Node* b @@
       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;
@@ @@ -1859,411 +1676,355 @@ void Mesh::DeleteLooseWalls(void) { @@
   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

@@ @@ -35,25 +35,28 @@ @@
 #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

@@ @@ -14,25 +14,28 @@ @@
  *  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

@@ @@ -152,67 +152,62 @@ void Node::DrawIndex(QGraphicsScene *c) @@
 		   ((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

@@ @@ -1859,24 +1859,25 @@ void Parameter::AssignValArrayToPar(cons @@
    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

@@ @@ -124,26 +124,27 @@ class Parameter { @@
   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

@@ @@ -782,24 +782,26 @@ void ParameterDialog::Reset(void) { @@
   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

@@ @@ -121,25 +121,28 @@ class ParameterDialog : public QDialog { @@
   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

@@ @@ -729,114 +729,108 @@ void Mesh::XMLSave(const char *docname, @@
   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);
@@ @@ -964,84 +958,78 @@ void Mesh::XMLReadGeometry(const xmlNode @@
+  }
   // 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");
@@ @@ -1087,51 +1075,49 @@ void Mesh::XMLReadNodes(xmlNode *root) @@
       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);
@@ @@ -1404,91 +1390,85 @@ void Mesh::XMLReadNodeSetsToNodes(xmlNod @@
       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);

0 comments (0 inline, 0 general)