/*
 *
 *  $Id$
 *
 *  This file is part of the Virtual Leaf.
 *
 *  The Virtual Leaf is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  The Virtual Leaf is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with the Virtual Leaf.  If not, see <http://www.gnu.org/licenses/>.
 *
 *  Copyright 2010 Roeland Merks.
 *
 */


#ifndef _NODESET_H_
#define _NODESET_H_

#include <algorithm>
#include <numeric>
#include <list>
#include <iterator>
#include "node.h"

class NodeSet : public list<Node *> {
  
 public:
  NodeSet(void) {
    done = false;
  }
  
  inline bool DoneP(void) { return done; }
  inline void ResetDone(void) { done = false; }
  
  void AddNode(Node * n) {
    push_back(n);
    n->node_set = this;
  }

  list <Cell *> getCells(void) {
    
    list<Cell *> cellset;
    
    for (list<Node *>::const_iterator i=begin();
	 i!=end();
	 i++) {
      transform ( (*i)->owners.begin(), (*i)->owners.end(), back_inserter( cellset ) , mem_fun_ref ( &Neighbor::getCell ));
    }
    
    cellset.sort();
    
    // remove all non-unique elements
    // (unique moves all unique elements to the front and returns an iterator to the end of the unique elements;
    // so we simply erase all remaining elements.
    cellset.erase(::unique(cellset.begin(), cellset.end() ), 
		  cellset.end() );
    
    // remove boundary_polygon
    cellset.erase( find_if ( cellset.begin(), cellset.end(), mem_fun( &Cell::BoundaryPolP  ) ) );
    return cellset;
  }

  void CleanUp(void) {
    
    // remove double Nodes from the set
    sort();
    erase(::unique(begin(), end() ), 
	      end() );
    
    
  }
  void print(ostream &os) const {
    transform( begin(), end(), ostream_iterator<int>( os, " " ), mem_fun( &Node::Index ) );
  }
  
  /*! Attempt a move over (rx, ry)
    reject if energetically unfavourable.
  */

  void AttemptMove(double rx, double ry) {
    
    done = true;
    // 1. Collect list of all attached to the nodes in the set
    list<Cell *> celllist = getCells();
    
    // 2. Sum the current energy of these cells
    double old_energy=0.;
    double sum_stiff = 0.;
    for ( list<Cell *>::const_iterator i = celllist.begin(); 
	  i!=celllist.end();
	  ++i ) {

      old_energy += (*i)->Energy();
      sum_stiff += (*i)->Stiffness();
    }
    
    // 3. (Temporarily) move the set's nodes.
    for ( list<Node *>::iterator n = begin();
	  n!=end();
	  ++n ) {
       
      (*n)->x+=rx;
      (*n)->y+=ry;
      // (*n)->z += rz;
       
    }
    
    // 4. Recalculate the energy
    double new_energy = 0.;
    for ( list<Cell *>::const_iterator i = celllist.begin(); 
	  i!=celllist.end();
	  ++i ) {
      
      new_energy += (*i)->Energy();
    }
    
    
    // 5. Accept or Reject DeltaH
	  double dh = new_energy - old_energy;
	  
	  list<int> new_areas;
	  
	 // cerr << "Nodeset says: dh = " << dh << " ...";
	  if (dh < -sum_stiff || RANDOM()<exp((-dh - sum_stiff)/par.T) ) {
		
		// ACCEPT
		// recalculate areas of cells
		for_each ( celllist.begin(), celllist.end(), mem_fun ( &Cell::RecalcArea ) );
		
		// recalculate integrals of cells
		for_each ( celllist.begin(), celllist.end(), mem_fun ( &Cell::SetIntegrals ) );
      
		
    } else {
      // REJECT
      
      // Move the set's nodes back to their original position
      for ( list<Node *>::iterator n = begin();
	    n!= end();
	    ++n ) {
	
	(*n)->x-=rx;
	(*n)->y-=ry;
	  }
      
      
    }
  }
  
  void XMLAdd(xmlNode *root) const;
  void XMLRead(xmlNode *root, Mesh *m);
 private:
  bool done;
};

ostream &operator<<(ostream &os, const NodeSet &ns);
#endif