#

#  $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.

#

CONFIG += qt

QMAKE_CXXFLAGS += -fexceptions

QMAKE_CXXFLAGS_DEBUG += -g3

QMAKE_CXXFLAGS_DEBUG += -DQDEBUG

#REACTIONS = reactions_auxin_growth.h 

#REACTIONS = reactions_meinhardt.h

#REACTIONS = reactions_pce_growth.h

DEFINES += QTGRAPHICS

DEFINES += REACTIONS_HEADER=$${REACTIONS}

DEFINES += REACTIONS_HEADER_STRING=\"$${REACTIONS}\"

DEFINES += FLEMING

PERLDIR = ./perl

BINDIR = ../bin

DESTDIR = $$BINDIR

TARGET = VirtualLeaf

TEMPLATE = app

PARTMPL = $${TARGET}par.tmpl

MAINSRC = $${TARGET}.cpp

QT -= network sql xml

QT += qt3support

#

# $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.

#

CONFIG += release

CONFIG -= debug

CONFIG += plugin

BINDIR = ../../bin

DEFINES = QTGRAPHICS # VLEAFPLUGIN

DESTDIR = $${BINDIR}/models

TARGET = test

HEADERS = ../simplugin.h $${TARGET}plugin.h  

QMAKE_CXXFLAGS += -fexceptions -I..

QMAKE_CXXFLAGS_DEBUG += -g3

QMAKE_CXXFLAGS_DEBUG += -DQDEBUG

QT += qt3support

SOURCES = $${TARGET}plugin.cpp

TEMPLATE = lib 

unix {

 LIBS += -Llib -lvleaf

 QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2

 QMAKE_LFLAGS += -fPIC

}

win32 {

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

	for (int i=0;i<NChem();i++) {

		daughter->chem[i]=chem[i];

	}

	//extern double auxin_account;

	//auxin_account += daughter->chem[0];

	for (int i=0;i<NChem();i++) {

		daughter->new_chem[i]=new_chem[i];

	}

	daughter->boundary=boundary;

	daughter->m=m;

	daughter->target_area=target_area/2.;

	target_area/=2;

	daughter->cellvec=cellvec;

//	daughter->BaseArea()  = base_area;

	// Division currently only works for convex cells: i.e. if the division line

	// intersects the cells at two points only.

	if (new_node_locations.size()!=2) {

		// Note: if you would add the possibility of dividing non-convex

		// cells, remember to update the code below. There are some

	new_chem=new double[NChem()];

	for (int i=0;i<NChem();i++) {

		new_chem[i]=0.;

	}

	boundary=None;

	index=(NCells()++);

	area=0.;

	target_area=1;

	target_length=0; //par.target_length;

	lambda_celllength = 0; //par.lambda_celllength;

	intgrl_xx=0.; intgrl_xy=0.; intgrl_yy=0.;

	intgrl_x=0.; intgrl_y=0.;

	source = false;

	source_conc = 0.;

	source_chem = 0;

	at_boundary=false;

	fixed = false;

	pin_fixed = false;

	stiffness = 0;

	marked = false;

	dead = false;

	div_counter=0;

	cell_type = 0;

	flag_for_divide = false;

}

CellBase::CellBase(double x,double y,double z) : QObject(), Vector(x,y,z)

{

#ifndef VLEAFPLUGIN

	if (static_data_members == 0) {

		static_data_members = new CellsStaticDatamembers();

	}

#endif

	chem=new double[NChem()];

	for (int i=0;i<NChem();i++) {

		chem[i]=0.;

	}

	new_chem=new double[NChem()];

	for (int i=0;i<NChem();i++) {

		new_chem[i]=0.;

	}

	boundary=None;

	area=0.;

	target_area=1;

	target_length=0; //par.target_length;

	lambda_celllength=0; // par.lambda_celllength;

	index=(NCells()++);

	intgrl_xx=0.; intgrl_xy=0.; intgrl_yy=0.;

	intgrl_x=0.; intgrl_y=0.;

	source = false;

	fixed = false;

	at_boundary=false;

	pin_fixed = false;

	stiffness = 0;

	marked=false;

	dead  = false;

	div_counter = 0;

	cell_type = 0;

	flag_for_divide = false;

}

CellBase::CellBase(const CellBase &src) :  Vector(src), QObject()

{

	chem=new double[NChem()];

	for (int i=0;i<NChem();i++) {

		chem[i]=src.chem[i];

	}

	new_chem=new double[NChem()];

	for (int i=0;i<NChem();i++) {

		new_chem[i]=src.new_chem[i];

	}

	boundary=src.boundary;

	area=src.area;

	target_length=src.target_length;

	lambda_celllength=src.lambda_celllength;

	intgrl_xx=src.intgrl_xx; intgrl_xy=src.intgrl_xy; intgrl_yy=src.intgrl_yy;

	intgrl_x=src.intgrl_x; intgrl_y=src.intgrl_y;

	target_area=src.target_area;

	index=src.index;

	nodes=src.nodes;

	neighbors=src.neighbors;

	walls=src.walls;

	source = src.source;

	fixed = src.fixed;

	source_conc = src.source_conc;

	source_chem = src.source_chem;

	cellvec = src.cellvec;

	at_boundary=src.at_boundary;

	pin_fixed = src.pin_fixed;

	stiffness = src.stiffness;

	marked = src.marked;

	dead = src.dead;

	cell_type = src.cell_type;

	div_counter = src.div_counter;

	flag_for_divide = src.flag_for_divide;

}

CellBase CellBase::operator=(const CellBase &src) {

	Vector::operator=(src);

	//  QObject::operator=(src);

	for (int i=0;i<NChem();i++) {

		chem[i]=src.chem[i];

	}

	for (int i=0;i<NChem();i++) {

		new_chem[i]=src.chem[i];

	}

	boundary=src.boundary;

	area=src.area;

	intgrl_xx=src.intgrl_xx; intgrl_xy=src.intgrl_xy; intgrl_yy=src.intgrl_yy;

	intgrl_x=src.intgrl_x; intgrl_y=src.intgrl_y;

	target_area=src.target_area;

	target_length=src.target_length;

	lambda_celllength=src.lambda_celllength;

	index=src.index;

	nodes=src.nodes;

	neighbors=src.neighbors;

	walls=src.walls;

	source = src.source;

	fixed = src.fixed;

	source_conc = src.source_conc;

	source_chem = src.source_chem;

	cellvec = src.cellvec;

	at_boundary=src.at_boundary;

	pin_fixed = src.pin_fixed;

	stiffness = src.stiffness;

	marked = src.marked;

	dead = src.dead;

	cell_type = src.cell_type;

	div_counter = src.div_counter;

	flag_for_divide = src.flag_for_divide;

	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;

}

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;

	double base_area;

};

class CellBase :  public QObject, public Vector 

{

	Q_OBJECT

	friend class Mesh;

	friend class CellInfo;

	friend class Node;

	friend class WallBase;

	friend class SimPluginInterface;

 public:

	CellBase(QObject *parent=0);

	CellBase(double x,double y,double z=0); // constructor

  virtual ~CellBase() {

	  delete[] chem;

	  delete[] new_chem;

	  //cerr << "CellBase " << index << " is dying. " << endl;

  }

	CellBase(const CellBase &src); // copy constructor

	virtual bool BoundaryPolP(void) const { return false; } 

	//  CellBase(const Vector &src); // not allowed (we cannot know to which mesh 

	/// the CellBase will belong...)

    CellBase operator=(const CellBase &src); // assignment operator

    CellBase operator=(const Vector &src);

    void SetChemical(int chem, double conc);

    inline void SetNewChem(int chem, double conc) { 

      new_chem[chem] = conc;

    }

    void SetSource(int chem, double conc) {

      source=true;

      source_chem = chem;

      source_conc = conc;

	}

    void UnfixNodes(void);

    void FixNodes(void);

      lambda_celllength = lambda_length;

    }

    inline double TargetArea(void) {

      return target_area;

    }

    inline void SetStiffness(double stiff) {

      stiffness = stiff;

    }

    inline double Stiffness(void) {

      return stiffness;

    }

    inline double EnlargeTargetArea(double da) {

      return target_area+=da;

    }

    inline double Area(void) const {

      return area;

    }

	inline void Divide(void) {

		flag_for_divide = true;

	}

    inline double Circumference(void) const {

      double sum=0.;

      for (list<Wall *>::const_iterator w=walls.begin();

	   w!=walls.end();

	   w++) {

	sum +=  (*w)->Length();

      }

      return sum;

    }

	QList<WallBase *> getWalls(void) {

		QList<WallBase *> wall_list;

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

			 i!=walls.end();

			 i++) {

			wall_list << *i;

		}

		return wall_list;

	}

  //  void XFigPrint(std::ostream &os) const;

    void Dump(ostream &os) const;

    // a (non-ordered) list of neighboring cells (actually I think the

    // introduction of ConstructWalls() has made these

    // lists ordered (clockwise), but I am not yet 100% sure...).

    list<CellBase *> neighbors;

    list<Wall *> walls;

    double *chem;

    double *new_chem;

    boundary_type boundary;

    mutable double area;

    double target_area;

    double target_length;

    double lambda_celllength;

    double stiffness; // stiffness like in Hogeweg (2000)

    bool fixed;

    bool pin_fixed;

    bool at_boundary; 

    bool dead; 

	bool flag_for_divide;

	int cell_type;

    //double length;

    //Vector meanflux;

    //int valence;

    // for length constraint

    mutable double intgrl_xx, intgrl_xy, intgrl_yy, intgrl_x, intgrl_y;

    bool source;

    Vector cellvec;

	// STATIC DATAMEMBERS MOVED TO CLASS

	static CellsStaticDatamembers *static_data_members;

	double source_conc;

    int source_chem;

    // PRIVATE MEMBER FUNCTIONS

    inline static void ClearNCells(void) {

      NCells()=0;

    }

    bool marked;

    int div_counter;

  //static double *derivs = 0; 

  // derivs is allocated by RungeKutta class.

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

	  plugin->CellDynamics(*c, &(derivs[i]));

	  i+=nchems;

  }

  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

	  //(*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 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) {

				(*i)->flag_for_divide=false;

			}

		}

	}

/*	template<class Op1, class Cont> void ExtractFromCells(Op1 f, Cont res) {

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

			 i!=cells.end();

			 i++) {

			*(res++) = ( f(*i) );

		}

	}*/

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