#

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

 matrix.h \

 mesh.h \

 miscq.h \

 modelcatalogue.h \

 Neighbor.h \

 node.h \

 nodeitem.h \

 nodeset.h \

 OptionFileDialog.h \

 output.h \

 parameter.h \

 pardialog.h \

 parse.h \

 pi.h \

 qcanvasarrow.h \

 random.h \

 rungekutta.h \

 simitembase.h \

 simplugin.h \

 sqr.h \

 tiny.h \

 transporterdialog.h \

 UniqueMessage.h \

 vector.h \

 wallbase.h \

 wall.h \

 wallitem.h \

 warning.h \

 xmlwrite.h \

 $${PARTMPL}

SOURCES += \

 apoplastitem.cpp \

 canvas.cpp \

 cellbase.cpp \

 cell.cpp \

 cellitem.cpp \

 forwardeuler.cpp \

 mainbase.cpp \

 matrix.cpp \

 mesh.cpp \

 miscq.cpp \

 modelcatalogue.cpp \

 Neighbor.cpp \

 node.cpp \

 nodeitem.cpp \

 nodeset.cpp \

 OptionFileDialog.cpp \

 output.cpp \

 parameter.cpp \

 pardialog.cpp \

 parse.cpp \

 random.cpp \

 rungekutta.cpp \

 simitembase.cpp \

 transporterdialog.cpp \

 UniqueMessage.cpp \

 vector.cpp \

 wallbase.cpp \

 wall.cpp \

 wallitem.cpp \

 warning.cpp \

 xmlwritecode.cpp \

 xmlwrite.cpp \

 $$MAINSRC

contains( TARGET, leaf_fleming ) {

 DEFINES += FLEMING	

}

contains(GRAPHICS, qwt) {

 #macx:LIBS += -L$$QWTDIR/lib -lqwt

 #win32:LIBS += -L$$QWTDIR/lib -lqwt5

 #LIBS += -L$$QWTDIR/lib -lqwt

 INCLUDEPATH += $$QWTDIR/include

 DEFINES += HAVE_QWT

 HEADERS += data_plot.h

 SOURCES += data_plot.cpp

}

contains( GRAPHICS, qt ) {

 message( "Building Qt executable" )

 QMAKE_CXXFLAGS += -DQTGRAPHICS # -fpermissive

}

contains( GRAPHICS, xfig ) {

 message("Building Xfig executable (background runnable).")

 QMAKE_CXXFLAGS += -DXFIGGRAPHICS

}

contains( GRAPHICS, x11 ) {

 !unix {

  error("X11 graphics only available on Unix systems.")

 }

 message("Building X11 executable")

 SOURCES += x11graph.cpp

 HEADERS += x11graph.h

 QMAKE_CXXFLAGS += -DX11GRAPHICS

 CONFIG -= qt

 CONFIG += x11

 unix:LIBS += -lpng

}

# finis

#

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

/*

 *

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

 *

 */

#include <string>

#include "cell.h"

#include "node.h"

#include "mesh.h"

#include "tiny.h"

#include "nodeset.h"

#include "cellitem.h"

#include "nodeitem.h"

#include "qcanvasarrow.h"

#include "parameter.h"

#include <QDebug>

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

extern Parameter par;

double Cell::factor=1.;

double Cell::offset[3]={0,0,0};

Cell::Cell(void) : CellBase() {

	m=0;

}

Cell::Cell(double x, double y, double z) : CellBase(x,y,z) {

	m=0;

}

Cell::Cell(const Cell &src) :  CellBase(src) {

	m=src.m;

}

bool Cell::Cmp(Cell *c) const { return this->Index() < c->Index(); }

bool Cell::Eq(Cell *c) const { return this->Index() == c->Index(); }

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

	CellBase::operator=(src);

	m=src.m;

	return *this;

}

//Cell(void) : CellBase() {}

void Cell::DivideOverAxis(Vector axis) {

	// Build a wall

	// ->  find the position of the wall

	// better look for intersection with a simple line intersection algorithm as below?

	// this leads to some exceptions: e.g. dividing a horizontal rectangle.

	// leaving it like this for the time being

	if (dead) return;

	Vector centroid=Centroid();

	double prev_cross_z=(axis * (centroid - *(nodes.back()) ) ).z ;

	ItList new_node_locations;

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

		 i!=nodes.end();

		 i++) {

		// cross product to detect position of division

		Vector cross = axis * (centroid - *(*i));

		if (cross.z * prev_cross_z < 0 ) {

			new_node_locations.push_back(i);

		} // else {

		//       //cerr << "cross.z * prev_cross_z = " << cross.z * prev_cross_z << endl;

		//     }

		prev_cross_z=cross.z;

	}

	DivideWalls(new_node_locations, centroid, centroid+axis);

}

double Cell::MeanArea(void) {

	return m->MeanArea();

}

void Cell::Apoptose(void) {

	// First kill walls

	cerr << "This is cell " << Index() << "\n";

	cerr << "Number of walls: " << walls.size() << endl;

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

		 w!=walls.end();

		 w++) {

		cerr << "Before apoptosis, wall " << (*w)->Index() << " says: c1 = " << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl;

	}

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

		 w!=walls.end();

		 w++) {

		bool illegal_flag = false;

		if ((*w)->c1 == (*w)->c2 ) 

			illegal_flag=true;

		if ((*w)->c1 == this) {

			// invert wall?

			(*w)->c1 = (*w)->c2;      

			(*w)->c2 = m->boundary_polygon;

			Node *n1 = (*w)->n1;

			(*w)->n1 = (*w)->n2;

			(*w)->n2 = n1;

		} else {

			(*w)->c2 = m->boundary_polygon;

		}

		if (illegal_flag && (*w)->c1==(*w)->c2) {

			cerr << "I created an illegal wall.\n";

		}

		if ( ((*w)->N1()->DeadP() || (*w)->N2()->DeadP()) ||

			((*w)->C1() == (*w)->C2() ) ){

			// kill wall

			cerr << "Killing wall.\n";

			(*w)->Kill();

			if ((*w)) {

				cerr << "Wall " << (*w)->Index() << " says: c1 = " << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl;

			}

			(*w)=0;

		} else {

			cerr << "Not killing wall.\n";

			cerr << "Wall " << (*w)->Index() << " says: c1 = " << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl;

		}

	}

	walls.remove(0);

	// Unregister me from my nodes, and delete the node if it no longer belongs to any cells

	list<Node *> superfluous_nodes;

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

		 n!=nodes.end();

		 n++) {

		Node &no(*(*n));

		// locate myself in the node's owner list

		list<Neighbor>::iterator cellpos;

		bool cell_found=false;

		for (list<Neighbor>::iterator nb=no.owners.begin();

			 nb!=no.owners.end();

			 nb++) {

			if (nb->cell == this) {

				cellpos = nb;

				cell_found = true;

				break;

			}

		}

		if (!cell_found) {

			// I think this cannot happen; but if I am wrong, unpredictable things would happen. So throw an exception.

			throw ("Cell not found in CellBase::Apoptose()\n\rPlease correct the code to handle this situation.");

		}

		Neighbor tmp = *cellpos;

		no.owners.erase(cellpos);

		// if node has no owners left, or only has a connection to special cell -1 (outside world), mark it as dead.

		if (no.owners.size()==0 || (no.owners.size()==1 && no.owners.front().cell->BoundaryPolP()) ) {

			no.MarkDead();

		} else {

			// register node with outside world

			if (find_if( no.owners.begin(), no.owners.end(), 

				     bind2nd ( mem_fun_ref(&Neighbor::CellEquals), m->boundary_polygon->Index() ) ) == no.owners.end() ) {

				tmp.cell = m->boundary_polygon;

				no.owners.push_back(tmp);

			}

		}

	}

	/*

	 // correct boundary polygon if this cell touches the boundary

	 // find the first living boundary node after a dead node

	 bool node_found = false;

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

	 n!=nodes.end();

	 n++) {

	 Node &no(*(*n));

	 if (no.DeadP()) {

	 list<Node *>::iterator first_node = n; 

	 if (++next_node == nodes.end()) first_node=nodes.begin();

	 if (!(*(*first_node)).DeadP() && ((*first_node)->boundary)) {

	 node_found=true;

	 break;

	 }

	 }

	 }

	 // locate it in the boundary_polygon

	 if (node_found) {

	 list<Node *>::iterator insert_it = find(mesh->boundary_polygon->nodes.begin(),

	 mesh->boundary_polygon->nodes.end(),

	 ++first_node);

	 if (insert_it!=owners.end()) {

	 if (insert_it==owners.end()) insert_it=owners.begin();

	 for (list<Node *>::iterator n=insert_it;

	 n!=nodes.end();

	 n++) {

	 Node &no(*(*n));

	 mesh->boundary_polygon->nodes.insert(

	 }

	 }

	 } */

	// mark cell as dead

	MarkDead();

}

void Cell::ConstructConnections(void) {

    // Tie up the nodes of this cell, assuming they are correctly ordered

    //cerr << "Constructing connections of cell " << index << endl;

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

		 i!=nodes.end();

		 i++) {

		//cerr << "Connecting node " << *i << endl;

		//cerr << "Node " << *i << endl << " = " << *(*i) << endl;

		// 1. Tidy up existing connections (which are part of this cell)

		if ((*i)->owners.size()>0) {

			list<Neighbor>::iterator neighb_with_this_cell=

			// remove myself from the neighbor list of the node

			find_if((*i)->owners.begin(),

					(*i)->owners.end(),

				bind2nd(mem_fun_ref( &Neighbor::CellEquals ),this->Index() )  );

			if (neighb_with_this_cell!=(*i)->owners.end()) 

				(*i)->owners.erase(neighb_with_this_cell);

		}

		Node *previous;

		if (i!=nodes.begin()) {

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

			previous_iterator--;

			previous=*previous_iterator;

		} else {

			previous=nodes.back();

		}

		Node *next;

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

		next_iterator++;

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

			next=nodes.front();

		} else {

			next=*next_iterator;

		}

		//cerr << "[" << *i << "]";

		//if (*i==10 || *i==11) {

		//cerr << "previous = " << previous << ", next = " << next << endl;

		//}

		//if (*i!=10 && *i!=11)

		//cerr << "Node " << *i << endl << " = " << *(*i) << endl;

		(*i)->owners.push_back( Neighbor( this, previous, next ) );

		// if (*i==50 || *i==51) {

		//cerr << "Node " << *i << ".size() = " << (*i)->owners.size() << endl;

		// }

    }

}

/*! \brief Divide the cell over the line v1-v2.

 \param v1: First vertex of line.

 \param v2: Second vertex of line.

 \param fixed_wall: If true: wall will be set to "fixed" (i.e. not motile)

 \return: true if the cell divided, false if not (i.e. no intersection between v1 and v2, and the cell)

 */

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