}

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