void operator() (Cell &c,std::ostream &os) const {

    os << "Cell " << c.index << " says: " << endl;

    os << "c.nodes.size() = " << c.nodes.size() << endl;

      cerr << (*i)->Index() << " ";

    }

    cerr << endl;

      NodeSet *node_set = new NodeSet;

	(*it)->setBrush(QBrush("purple"));

	Cell &c=(*it)->getCell();

	// sometimes the cell hasn't properly divided yet before the

	mesh.CutAwayBelowLine( startpoint, endpoint ); 

	// Correct flags of nodeset

	  (*i)->SetSAM();

	  (*i)->SetBoundary();

	}

	//   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();

	  (*c)->SetBoundary(Cell::SAM);

	}

  ItList new_node_locations;

    // cross product to detect position of division

    Vector cross = axis * (centroid - *(*i));

	     << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl;

#endif

  }

    bool illegal_flag = false;

    if ((*w)->c1 == (*w)->c2 ) 

  // Unregister me from my nodes, and delete the node if it no longer belongs to any cells

  list<Node *> superfluous_nodes;

    Node &no(*(*n));

    // locate myself in the node's owner list

    list<Neighbor>::iterator cellpos;

    bool cell_found=false;

      if (nb->cell == this) {

	cellpos = nb;

	cell_found = true;

  //cerr << "Constructing connections of cell " << index << endl;

    //cerr << "Connecting node " << *i << endl;

    //cerr << "Node " << *i << endl << " = " << *(*i) << endl;

#ifdef QDEBUG

  qDebug() << "Cell " << Index() << " is doing DivideOverGivenLine" << endl;

#endif

    Vector v3 = *(*i);

    list<Node *>::iterator nb=i;

      new_node_locations.push_back(nb);

    } 

#ifdef QDEBUG

  if (new_node_locations.size()<2) { 

    return false;

  }

#endif

}

// Core division procedure

  // construct new Nodes at the intersection points

  // unless they coincide with existing points

    // intersection between division axis

    // and line from this node to its predecessor

      owners.sort( mem_fun_ref( &Neighbor::Cmp ) );

      list<Neighbor>::iterator c;

	c=adjacent_find(c,owners.end(),neighbor_cell_eq);

	if (c->cell->Index() != this->Index() || c==owners.end()) break;

      }

  // Move Walls to daughter cell

  list <Wall *> copy_walls = walls;

    //cerr << "Doing wall, before:  " << **w << endl;

  // now reconstruct neighbor list for all "broken" neighbors

    ((Cell *)(*i))->ConstructNeighborList();

  }

// Move the whole cell

void Cell::Move(const Vector T) {

    *(*i)+=T;

  }

}

  // the length of these edges will change

  double old_length=0.;

      if (n->getCell()!=this) {

	length_edges.push_back( pair <Node *,Node *> (*i, n->nb1) );

	length_edges.push_back( pair <Node *,Node *> (*i, n->nb2) );

  // calculate area energy difference of neighboring cells

  // (this cells' shape remains unchanged)

  double old_area_energy=0., old_length_energy=0.;

    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.;

    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;

  }

    // update areas of cells

    //cerr << "neighbors: ";

    list<CellBase *>::const_iterator nb_it = neighbors.begin();

      ((Cell *)(*nb_it))->area = *ar_it;

      (*nb_it)->SetIntegrals(); 

    }

  double energy = 0.;

  double length_contribution = 0.;

      if (n->getCell()==this) {

	length_contribution += 

	  DSQR(Node::target_length-(*(*i)-*(n->nb2)).Norm());

      }

    }

  }

  // method used: http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/

    list<Node *>::const_iterator j=i; 

    ++j;

      {

	Vector v1 = *(*i);

  neighbor_of_moving_node[1]=*( *nb ); 

    for (int j=0;j<2;j++) { // loop over the two neighbors of moving node

      list<Node *>::const_iterator nb=i;

      nb++;

  // Get "corner points; i.e. nodes where more than 2 cells are connected

  list<Node *> corner_points;

    // look for nodes belonging to >2 cells

    if ((*i)->owners.size()>2) {

  list<Node *>::const_iterator nb = (--corner_points.end());

  // loop over list, 

    if (nb==corner_points.end()) nb=corner_points.begin();

    // add owning cells to a list

    list<Cell *> owning_cells;

    Node &n(*(*i));

      owning_cells.push_back(j->cell);

    }

    Node &n2(*(*nb));

      owning_cells.push_back(j->cell);

    }

    // find duplicates

    vector<Cell *> duplicates;

    list<Cell *>::const_iterator prevj = (--owning_cells.end());

    }

    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()) 

  QPolygonF pa(nodes.size());

  int cc=0;

    Node *i=*n;

  }

  p->setPolygon(pa);

{

  // loop over cell edges

    // leaf cannot take up chemicals from environment ("no flux boundary")

    if ((*i)->c2->BoundaryPolP()) continue;

  QPolygonF pa(nodes.size());

  int cc=0;

    Node *i=*n;

    pa[cc++] = QPoint((int)((offset[0]+i->x)*factor),

void Cell::DrawNodes(QGraphicsScene *c) const {

    Node *i=*n;

    NodeItem *item = new NodeItem ( &(*i), c );

    item->setColor();

    item->setZValue(5);

    item->show();

  }

}

void Cell::SetWallLengths(void)

{

    // Step 1: find the path of nodes leading along the Wall.

    // A Wall often represents a curved cell wall: we want the total

    double sum_length = 0.;

    // Now, walk to the second node of the edge in the list of nodes

      list<Node *>::const_iterator prev_n = n; 

      --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(); 

  walls.push_back( w );

  // Add wall to Mesh's list if it isn't there yet

  }

}

{

  // remove wall from Mesh's list

  // 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);

}

    error << "SetChemical: value ch = " << ch << " is out of range\n";

    throw error.str().c_str();

  }

    (*w)->setTransporter(this, ch, conc);

  }

}

  os << " [ area = " << area << " ]";

  os << " [ walls = ";

    os << (*i)->n1->Index() << " -> " << (*i)->n2->Index() << ", " <<  (*i)->c1->Index() << " | " << (*i)->c2->Index() << ", ";

  }

  os << " ] ";

  double loc_area=0.;

    list<Node *>::const_iterator i_plus_1=i; i_plus_1++;

    if (i_plus_1==nodes.end())

  double area=0.;

  double integral_x_dxdy=0.,integral_y_dxdy=0.;

    list<Node *>::const_iterator i_plus_1=i; i_plus_1++;

    if (i_plus_1==nodes.end())

  list<Node *>::const_iterator nb;

  list<Node *>::const_iterator i=nodes.begin();

    nb = i; nb++; if (nb==nodes.end()) nb=nodes.begin();

  list<Node *>::const_iterator nb;

  list<Node *>::const_iterator i=nodes.begin();

    nb = i; nb++; if (nb==nodes.end()) nb=nodes.begin();

  os << " " << source_conc << " " << source_chem;

  os << endl;

}

void CellBase::UnfixNodes(void)

{

    (*i)->Unfix();

  }

}

void CellBase::FixNodes(void)

{

    (*i)->Fix();

  }

}

// returns true if cell is at border

void Mesh::PerturbChem(int chemnum, double range) {

    (*i)->chem[chemnum] += range*(RANDOM()-0.5);

    if ((*i)->chem[chemnum]<0.) (*i)->chem[chemnum]=0.;

  }

    IncreaseCellCapacityIfNecessary();

    vector <Cell *> current_cells = cells;

      (*j)->DivideOverAxis(axis);

    }

    axis=axis.Perp2D();

  axis=axis.Perp2D();

  vector <Cell *> current_cells = cells;

    (*j)->DivideOverAxis(axis);

  }

  double sum_l=0; int n_l=0;

    list<Node *>::const_iterator nb=i; nb++; 

    if (nb==cell->nodes.end()) 

      nb=cell->nodes.begin();

  // a bit more tension

  Node::target_length/=2;

  c->SetIntegrals(); 

  c->at_boundary=true;

  return *c;

  // get position of the (n/4)'th and (3*(n/4))'th node.

  list<Node *>::reverse_iterator it_n1=circle->nodes.rbegin();

    it_n1++;

  it_n1--;

  list<Node *>::reverse_iterator it_n2=--circle->nodes.rend();

  Cell *petiole=AddCell(new Cell());

  list<Node *>::reverse_iterator i=it_n1; i++;

		  nodes[(*it_n2)->Index() + (((*i)->Index()-(*it_n2)->Index()) + 1)%nnodes],

		  nodes[(*it_n2)->Index() + (((*i)->Index()-(*it_n2)->Index())-1+nnodes)%nnodes]);

  }

  (*it_n2)->boundary=true;

#ifdef QDEBUG  

  qDebug() << petiole << endl;

#endif

  (*it_n1)->boundary=true;

    }

  }

  // make petiole solid

    (*i)->Fix();

  }

  petiole->Fix();

  LowerLeft = **nodes.begin();

  UpperRight = **nodes.begin();

    if ((*c)->x < LowerLeft.x)

      LowerLeft.x = (*c)->x;

    if ((*c)->y < LowerLeft.y)

void Mesh::Clear(void) {

  // clear nodes

    delete *i;

  }

  node_insertion_queue.clear();

  // Clear NodeSets

    delete *i;

  }

  // clear cells

    delete *i;

  }

  delete boundary_polygon;

  // Clear walls