// relocated including the current Cell (i.e. the value of *this)

    // calling "Mesh::IncreaseCapacityIfNecessary" (from another

    // object than Cell, e.g. Mesh) before entering Divide will solve

    // this issue (solved now).

    cells.reserve(2);

    nodes.reserve(500);

    time = 0.;

    plugin = 0;

    boundary_polygon=0;

  };

  ~Mesh(void) {

    if (boundary_polygon) {

      delete boundary_polygon;

      boundary_polygon=0;

    }

  };

  void Clean(void);

  Cell &EllipticCell(double xc, double yc, double ra, double rb, int nnodes=10, double rotation=0);

  Cell &CircularCell(double xc, double yc, double r, int nnodes=10) {

    return EllipticCell(xc, yc, r, r, nnodes, 0);

  }

  Cell &LeafPrimordium(int n, double pet_length);

  Cell &LeafPrimordium2(int n);

  Cell *RectangularCell(const Vector ll, const Vector ur, double rotation = 0);

  void CellFiles(const Vector ll, const Vector ur);

  inline Cell &getCell(int i) {

    if ((unsigned)i<cells.size())

      return *cells[i];

    else {

#ifdef QDEBUG

      qDebug() << i << endl;

      qDebug() << "size is " << cells.size() << endl;

#endif

      abort();

    }

  }

  inline Node &getNode(int i) {

    return *nodes[i];    

  }

  //double Diffusion(void);

  inline int size(void) {

    return cells.size();

  }

  inline int nnodes(void) {

    return nodes.size();

  }

  template<class Op> void LoopCells(Op f) {

    for (vector <Cell *>::iterator i=cells.begin();

	 i!=cells.end();

	 i++) {

      f(**i);

    }

  }

  template<class Op> void LoopWalls(Op f) {

    for (list <Wall *>::iterator i=walls.begin();

	 i!=walls.end();

	 i++) {

      f(**i);

    }

  }

  // if the amount of cells might increase, during looping, use this template

  template<class Op> void LoopCurrentCells(Op f) {

    vector<Cell *> current_cells = cells;

    for (vector <Cell *>::iterator i=current_cells.begin();

	 i!=current_cells.end();

	 i++) {

      f(**i);

    }

  }

  template<class Op> void LoopNodes(Op f) {

    for (vector<Node *>::iterator i=nodes.begin();

	 i!=nodes.end();

	 i++) {

      f(**i); 

    }

  }

  template<class Op> void RandomlyLoopNodes(Op f) {

    MyUrand r(shuffled_nodes.size());

    random_shuffle(shuffled_nodes.begin(),shuffled_nodes.end(),r);

    for (vector<Node *>::const_iterator i=shuffled_nodes.begin();

	 i!=shuffled_nodes.end();

	 i++) {

    }

  }

  template<class Op> void RandomlyLoopCells(Op f) {

    MyUrand r(shuffled_cells.size());

    random_shuffle(shuffled_cells.begin(),shuffled_cells.end(),r);

    for (vector<Cell *>::const_iterator i=shuffled_cells.begin();

	 i!=shuffled_cells.end();

	 i++) {

    }

  }

  template<class Op1, class Op2> void LoopCells(Op1 f, Op2 &g) {

    for (vector<Cell *>::iterator i=cells.begin();

	 i!=cells.end();

	 i++) {

      f(**i,g); 

    }

  }

  template<class Op1, class Op2, class Op3> void LoopCells(Op1 f, Op2 &g, Op3 &h) {

    for (vector<Cell *>::iterator i=cells.begin();

	 i!=cells.end();

	 i++) {

      f(**i,g,h); 

    }

  }

  void DoCellHouseKeeping(void) {

    vector<Cell *> current_cells = cells;

    for (vector<Cell *>::iterator i = current_cells.begin();

	 i != current_cells.end();

	 i ++) {

      plugin->CellHouseKeeping(*i);

      // Call functions of Cell that cannot be called from CellBase, including Division

      if ((*i)->flag_for_divide) {

	if ((*i)->division_axis) {

	  (*i)->DivideOverAxis(*(*i)->division_axis);

	  delete (*i)->division_axis;

	  (*i)->division_axis = 0;

	} else {

	  (*i)->Divide();

	}

	(*i)->flag_for_divide=false;

      }

    }

  }

  // Apply "f" to cell i

  // i.e. this is an adapter which allows you to call a function

  // operating on Cell on its numeric index index

  template<class Op> void cell_index_adapter(Op f,int i) {

    f(cells[i]);

  }

  double DisplaceNodes(void);

  void BoundingBox(Vector &LowerLeft, Vector &UpperRight);

  int NEqs(void) {     int nwalls = walls.size();

    int ncells =cells.size();

    int nchems = Cell::NChem();

    // two eqs per chemical for each walls, and one eq per chemical for each cell

    // This is for generality. For a specific model you may optimize

    // this by removing superfluous (empty) equations.

    int neqs = 2 * nwalls * nchems + ncells * nchems;

    return neqs;

  }

  void IncreaseCellCapacityIfNecessary(void) {

    return;

    // cerr << "Entering Mesh::IncreaseCellCapacityIfNecessary \n";

    // make sure we always have enough space 

    // to have each cell divide at least once

    //

    // Note that we must do this, because Cell::Divide pushes a new Cell

    // onto Mesh::cells. As a result, Mesh::cells might be relocated 

    // if we are _within_ a Cell object: i.e. pointer "this" will be changed!!

    // 

    // An alternative solution could be to make "Mesh::cells" a list,

    // but this won't work because we need random access for 

    // the Monte Carlo algorithm.

    if (2*cells.size()>cells.capacity()) {

      cerr << "Increasing capacity to "  << 2*cells.capacity() << endl;

      cerr << "Current capacity is " << cells.capacity() << endl;

      cells.reserve(cells.capacity()*2);

    }

  }

  void ReserveMoreCells(int n) {

    if (nodes.size()+n>nodes.capacity()) {

      nodes.reserve(size()+n);

    }

  }

  double Area(void);

  double MeanArea(void) {

    double sum=0.;

    for (vector<Cell *>::const_iterator i=cells.begin();

	 i!=cells.end();

	 i++) {

      sum+=(*i)->Area();

    }