#

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

!win32 {

 GRAPHICS = qt #qwt

}

win32 {

 CONFIG += console

 LIBXML2DIR = C:\libxml2

 LIBICONVDIR = C:\libiconv

 LIBZDIR = C:\libz

 system(DEL parameter.cpp parameter.h) 

 GRAPHICS = qt 

 RC_FILE = VirtualLeaf.rc

 QMAKE_CXXFLAGS += -DLIBXML_STATIC

 QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include

 QMAKE_POST_LINK = "\

  C:\Bin\cp release\VirtualLeaf.exe \

  C:\Qt\4.5.3\bin\Qt3Support4.dll \

  C:\Qt\4.5.3\bin\QtGui4.dll \

  C:\Qt\4.5.3\bin\QtSql4.dll \

  C:\Qt\4.5.3\bin\QtXml4.dll \

  C:\Qt\4.5.3\bin\QtCore4.dll \

  C:\Qt\4.5.3\bin\QtNetwork4.dll \

  C:\Qt\4.5.3\bin\QtSvg4.dll \

  C:\bin\iconv.dll \

  C:\bin\libxml2.dll \

  C:\bin\zlib1.dll \

  C:\MinGW\bin\mingwm10.dll \

  $${DESTDIR}"

 LIBS += -L$${LIBXML2DIR}\lib -lxml2 -L$${LIBICONVDIR}\lib -L$${LIBZDIR}\lib  -lz -lm -lwsock32 -liconv

}

# Application icons

macx {

 ICON = leaficon.icns

 # make sure that the executable can find the libqwt share library

 QMAKE_POST_LINK = "mkdir $${DESTDIR}/$${TARGET}.app/Contents/Frameworks; \

  cp /usr/local/qwt/lib/libqwt.dylib $${DESTDIR}/$${TARGET}.app/Contents/Frameworks/.; \

  #install_name_tool -change libqwt.5.dylib $$QWTDIR/lib/libqwt.dylib $${DESTDIR}/$${TARGET}.app/Contents/MacOS/$${TARGET}; \

  install_name_tool -id @executable_path/../Frameworks/libqwt.dylib $${DESTDIR}/$${TARGET}.app/Contents/Frameworks/libqwt.dylib; \

  install_name_tool -change libqwt.5.dylib @executable_path/../Frameworks/libqwt.dylib $${DESTDIR}/$${TARGET}.app/Contents/MacOS/$${TARGET};\

  cp leaficon.icns $${DESTDIR}/$${TARGET}.app; "

}

macx:release {

 LIBS+= -dead_strip

}

unix {

 system(rm -f parameter.cpp parameter.h) 

 CC = /usr/bin/gcc 

 QWTDIR = /ufs/guravage/opt/qwt-5.2.1-svn

 QMAKE_LIBDIR += $$QWTDIR/lib 

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

 QMAKE_LFLAGS += -fPIC

 LIBS += -lxml2 -lz -lm 

}

system(perl $$PERLDIR/make_parameter_source.pl $$PARTMPL)

system(perl $$PERLDIR/make_pardialog_source.pl $$PARTMPL)

#system(perl $$PERLDIR/make_xmlwritecode.pl -h $$REACTIONS)

# Input

HEADERS += \

 apoplastitem.h \

 canvas.h \

 cellbase.h \

 cell.h \

 cellitem.h \

 forwardeuler.h \

 infobar.h \

 mainbase.h \

 mainbase.h \

#

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

 LIBXML2DIR = C:\libxml2

 LIBICONVDIR = C:\libiconv

 LIBZDIR = C:\libz

 LIBS += -Llib -lvleaf

 QMAKE_CXXFLAGS += -DLIBXML_STATIC

 QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include

}

# finis

bool Cell::DivideOverGivenLine(const Vector v1, const Vector v2, bool fix_cellwall, NodeSet *node_set ) {

	if (dead) return false;

	// check each edge for intersection with the line

	ItList new_node_locations;

	cerr << "Cell " << Index() << " is doing DivideOverGivenLine \n";

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

		 i!=nodes.end();

		 i++) {

		Vector v3 = *(*i);

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

		nb++;

		if (nb == nodes.end()) {

			nb = nodes.begin();

		}

		Vector v4 = *(*nb);

		double denominator = 

		(v4.y - v3.y)*(v2.x - v1.x) - (v4.x - v3.x)*(v2.y - v1.y);

		double ua = 

		((v4.x - v3.x)*(v1.y - v3.y) - (v4.y - v3.y)*(v1.x -v3.x))/denominator;

		double ub = 

		((v2.x - v1.x)*(v1.y-v3.y) - (v2.y- v1.y)*(v1.x - v3.x))/denominator;

		/* double intersec_x = v1.x + ua*(v2.x-v1.x);

		 double intersec_y = v1.y + ua*(v2.y-v1.y);*/

		//cerr << "Edge " << *i << " to " << *nb << ": ua = " << ua << ", ub = " << ub << ":  ";

		// this construction with "TINY" should simulate open/closed interval <0,1]

		if ( ( TINY < ua && ua < 1.+TINY ) && ( TINY < ub && ub < 1.+TINY ) ) {

			// yes, intersection detected. Push the location to the list of iterators

			new_node_locations.push_back(nb);

		} 

	}

	if (new_node_locations.size()<2) { 

		cerr << "Line does not intersect with two edges of Cell " << Index() << endl;

		cerr << "new_node_locations.size() = " << new_node_locations.size() << endl;

		return false;

	}

	ItList::iterator i = new_node_locations.begin();

	list< Node *>::iterator j;

	cerr << "-------------------------------\n";

	cerr << "Location of new nodes: " << (**i)->Index() << " and ";

	++i;

	j = *i; 

	if (j==nodes.begin()) j=nodes.end(); j--;

	cerr << (*j)->Index() << endl;

	cerr << "-------------------------------\n";

	if ( **new_node_locations.begin() == *j ) {

		cerr << "Rejecting proposed division (cutting off zero area).\n";

		return false;

	}

	DivideWalls(new_node_locations, v1, v2, fix_cellwall, node_set);

	return true;

}

// Core division procedure

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

		// fixed-size arrays over there!

		cerr << "Warning in Cell::Division: division of non-convex cells not fully implemented" << endl;

		// Reject the daughter cell and decrement the amount of cells

		// again. We can do this here because it is the last cell added.

		// Never, ever try to fully erase a cell elsewhere, because we

		// make heavy use of cell indices in this project; if you erase a

		// Cell somewhere in the middle of Mesh::Cells the indices will

		// get totally messed up...! (e.g. the indices used in Nodes::cells)

		cerr << "new_node_locations.size() = " << new_node_locations.size() <<endl;

		cerr << "daughter->index = " << daughter->index << endl;

		cerr << "cells.size() = " << m->cells.size() << endl;

		m->cells.pop_back();

		Cell::NCells()--;

		m->shuffled_cells.pop_back();

		return;

	}

	// We can be sure we only need two positions here because divisions

	// of non-convex cells are rejected above.

	Vector new_node[2];

	Node *new_node_ind[2];

	int new_node_flag[2];

	Edge div_edges[2];

	int nnc=0;

	Wall *div_wall[4];

	double orig_length[2];

	for (int i=0;i<4;i++) { div_wall[i]=0; orig_length[i/2] = 0.; }

	// construct new Nodes at the intersection points

	// unless they coincide with existing points

	for ( ItList::const_iterator i=new_node_locations.begin();

		 i!=new_node_locations.end();

		 i++) {

		// intersection between division axis

		// and line from this node to its predecessor

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

		Vector v1 = from;

		Vector v2 = to;

		Vector v3 = *(**i);

		// get previous node

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

		if (nb == nodes.begin()) {

			nb = nodes.end();

		} 

		nb--;

		Vector v4=*( *nb ); 

		double denominator = 

		(v4.y - v3.y)*(v2.x - v1.x) - (v4.x - v3.x)*(v2.y - v1.y);

		double ua = 

		((v4.x - v3.x)*(v1.y - v3.y) - (v4.y - v3.y)*(v1.x -v3.x))/denominator;

		double intersec_x = v1.x + ua*(v2.x-v1.x);

		double intersec_y = v1.y + ua*(v2.y-v1.y);

		// construct a new node at intersec

		// we construct a vector temporarily,

		// until we are sure we are going to keep the node...

		// Node count is damaged if we construct superfluous nodes

		Vector *n=new Vector(intersec_x,intersec_y,0);

 *  (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.

 *

 */

#include <cmath>

#include <string>

#include <sstream>

#include <vector>

#include <algorithm>

#include <functional>

#ifdef QTGRAPHICS

#include <QGraphicsScene>

#include <qpainter.h>

#include <qcolor.h>

#include <qfont.h>

#include <qwidget.h>

//Added by qt3to4:

#include <Q3PointArray>

#include <fstream>

#include "nodeitem.h"

#include "cellitem.h"

#include "qcanvasarrow.h"

#endif

#include "nodeset.h"

//#include "cond_operator.h"

#include "cellbase.h"

//#include "node.h"

#include "wall.h"

#include "random.h"

#include "parameter.h" 

#include "mesh.h"

#include "sqr.h"

#include "tiny.h"

static const std::string _module_id("$Id$");

extern Parameter par;

const char* CellBase::boundary_type_names[4] = {"None", "NoFlux", "SourceSink", "SAM"};

// These statics have moved to class "CellsStaticDatamembers"

//double CellBase::static_base_area = 0.;

//int CellBase::ncells=0;

//int CellBase::NChem()=0;

#ifndef VLEAFPLUGIN

CellsStaticDatamembers *CellBase::static_data_members = new CellsStaticDatamembers();

#else

CellsStaticDatamembers *CellBase::static_data_members = 0;

#endif

CellBase::CellBase(QObject *parent) : 

QObject(parent),

Vector()

{

	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;

	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;

	for (list<Node *>::const_iterator i =  nodes.begin(); i!=nodes.end(); i++) {

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

	}

	os << " } ";

	for (list<Wall *>::const_iterator i =  walls.begin(); i!=walls.end(); i++) {

		(*i)->print(os);

		os << ", ";

	} 

	os << endl;

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

	os << " [ walls = ";

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

		 i!=walls.end();

		 i++) {

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

	}

	os << " ] ";

	os << "div_counter = " << div_counter << endl;

	os << "cell_type = " << cell_type << endl;

	os << endl;

	return os;

}

ostream &operator<<(ostream &os, const CellBase &c) {

	c.print(os);

	return os;

}

double CellBase::CalcArea(void) const {

	double loc_area=0.;

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

		 i!=(nodes.end());

		 i++) {

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

		if (i_plus_1==nodes.end())

			i_plus_1=nodes.begin();

		loc_area+= (*i)->x * (*i_plus_1)->y;

		loc_area-= (*i_plus_1)->x * (*i)->y;

	}

	// http://technology.niagarac.on.ca/courses/ctec1335/docs/arrays2.pdf	

	//return loc_area/2.0; 

	return fabs(loc_area)/2.0; 

} 

Vector CellBase::Centroid(void) const {

	double area=0.;

	double integral_x_dxdy=0.,integral_y_dxdy=0.;

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

		 i!=(nodes.end());

		 i++) {

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

		if (i_plus_1==nodes.end())

			i_plus_1=nodes.begin();

		area+= (*i)->x * (*i_plus_1)->y;

		area-= (*i_plus_1)->x * (*i)->y;

		integral_x_dxdy+=

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

 *

 */

// CellBase derives from Vector, where Vector is simply used as a Vertex

#ifndef _CELLBASE_H_

#define _CELLBASE_H_

#include <list>

#include <vector>

#include <iostream>

#include <QString>

#include <QDebug>

#include "vector.h"

#include "parameter.h"

#include "wall.h"

#include "warning.h"

#include "assert.h"

extern Parameter par;

using namespace std;

class Mesh;

class Node;

class CellBase;

class NodeSet;

struct ParentInfo {

	Vector polarization;

	double PINmembrane;

	double PINendosome;

};

// We need a little trick here, to make sure the plugin and the main application will see the same static datamembers

// otherwise each have their own instantation.

// My solution is as follow. I collect all original statics in a class. The main application instantiates it and

// has a static pointer to it. After loading the plugin I set a static pointer to the same class 

class CellsStaticDatamembers {

public:

	CellsStaticDatamembers(void) {

		ncells = 0;

		nchem = 0;

		base_area = 0.;

		cerr << "Constructor of CellsStaticDatamembers\n";

	}

	~CellsStaticDatamembers() {

		cerr << "Oops! Desctructor of CellsStaticDatamembers called\n";

	}

	int ncells;

	int nchem;

	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;