Changeset - c8fbdf1e4153
[Not reviewed]
default
0 7 0
Michael Guravage - 15 years ago 2010-05-28 13:53:39
michael.guravage@cwi.nl
Roeland's recent improvements - especially a more robust Mesh::Derivatives().

--
user: Michael Guravage <michael.guravage@cwi.nl>
branch 'default'
changed src/VirtualLeaf.pro
changed src/build_models/plugin_test.pro
changed src/cell.cpp
changed src/cellbase.cpp
changed src/cellbase.h
changed src/mesh.cpp
changed src/mesh.h
7 files changed with 36 insertions and 6 deletions:
0 comments (0 inline, 0 general)
src/VirtualLeaf.pro
Show inline comments
 
#
 
#  $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 += release
 
CONFIG -= debug
 
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
src/build_models/plugin_test.pro
Show inline comments
 
#
 
# $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
src/cell.cpp
Show inline comments
 
@@ -226,384 +226,385 @@ void Cell::Apoptose(void) {
 
	 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];
 
	}
 
	
 
	daughter->cell_type = cell_type;
 
	//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);
 
		
 
		div_edges[nnc].first=*nb;
 
		div_edges[nnc].second=**i;
 
		
 
		// Insert this new Node if it is far enough (5% of element length)
 
		// from one of the two existing nodes, else use existing node
 
		//
 
		// old, fixed value was: par.collapse_node_threshold = 0.05
 
		double collapse_node_threshold = 0.05;
 
#ifdef FLEMING
 
		collapse_node_threshold = par.collapse_node_threshold;
 
#endif
 
		
 
		double elem_length = ( (*(**i)) - (*(*nb)) ).Norm();
 
		if ( ( *(**i) - *n ).Norm() < collapse_node_threshold  * elem_length ) {
 
			new_node_flag[nnc]=1;
 
			new_node[nnc] = *(**i);
 
			new_node_ind[nnc] = **i;
 
			//cerr << **i << "\n" ;
 
		} else 
 
			if ( (*(*nb) - *n).Norm() < collapse_node_threshold * elem_length ) {
 
				new_node_flag[nnc]=2;
 
				new_node[nnc] = *(*nb);
 
				new_node_ind[nnc] = *nb;
 
			} else {
 
				new_node_flag[nnc]=0;
 
				new_node[nnc] = *n;
 
			}
 
		
 
		nnc++;
 
		delete n;
 
	}
 
	
 
	
 
	for (int i=0;i<2;i++) {
 
		
 
		Cell *neighbor_cell=0; // we need this to split up the "Wall" objects.
 
		
 
		// for both divided edges: 
 
		//      insert its new node into all cells that own the divided edge
 
		// but only if it really is a new node:
 
		if (new_node_flag[i]!=0) {
 
			if (fix_cellwall) {
 
				(new_node_ind[i])->fixed = true;
 
				
 
				// all this we'll do later for the node set :-)
 
				/* (new_node_ind[i])->boundary = true;
 
				 (new_node_ind[i])->sam = true;
 
				 boundary = SAM;
 
				 daughter->boundary = SAM;
 
				 boundary_touched_flag = true;
 
				 */ 
 
			}
 
			
 
		} else {
 
			
 
			// (Construct a list of all owners:)
 
			// really construct the new node (if this is a new node)
 
			new_node_ind[i] = 
 
			m->AddNode(new Node (new_node[i]) );
 
			
 
			
 
			
 
			// if a new node is inserted into a fixed edge (i.e. in the petiole)
 
			// make the new node fixed as well
 
			(new_node_ind[i])->fixed = (div_edges[i].first)->fixed &&
 
			(div_edges[i].second)->fixed;
 
			
 
			// Insert Node into NodeSet if the div_edge is part of it.
 
			if (
 
				(div_edges[i].first->node_set && div_edges[i].second->node_set) &&
 
				(div_edges[i].first->node_set == div_edges[i].second->node_set))
 
			{
 
				//cerr << "Inserting node into node set\n";
 
				div_edges[i].first->node_set->AddNode( new_node_ind[i] );
 
			}
 
			
 
			// if the new wall should be fixed (i.e. immobile, or moving as
 
			// solid body), make it so, and make it part of the boundary. Using
 
			// this to make a nice initial condition by cutting off part of a
 
			// growing leaf.
 
			
 
			if (fix_cellwall) {
 
				(new_node_ind[i])->fixed = true;
 
				
 
				// All this we'll do later for the node set only
 
				/* (new_node_ind[i])->boundary = true;
 
				 (new_node_ind[i])->sam = true;
 
				 boundary_touched_flag = true;
 
				 boundary = SAM;
 
				 daughter->boundary = SAM;*/
 
			}
 
			
 
			// if new node is inserted into the boundary
 
			// it will be part of the boundary, too
 

	
 
			new_node_ind[i]->UnsetBoundary();
src/cellbase.cpp
Show inline comments
 
/*
 
 *
 
 *  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 <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;
 
	division_axis = 0;
 
}
 

	
 

	
 
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;
 
	division_axis = 0;
 

	
 
}
 

	
 
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;
 
	</i>
 
	division_axis = src.division_axis;
 
}
 

	
 

	
 
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;
 
	division_axis = src.division_axis;
 
	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+=
 
			((*i_plus_1)->x+(*i)->x)*
 
			((*i)->x*(*i_plus_1)->y-
 
			 (*i_plus_1)->x*(*i)->y);
 
		integral_y_dxdy+=
 
			((*i_plus_1)->y+(*i)->y)*
 
			((*i)->x*(*i_plus_1)->y-
 
			 (*i_plus_1)->x*(*i)->y);
 
	}
 
    
 
	//area/=2.0;
 
	area = fabs(area)/2.0;
 
	
 
	integral_x_dxdy/=6.;
 
	integral_y_dxdy/=6.;
 
	
 
	Vector centroid(integral_x_dxdy,integral_y_dxdy,0);
 
	centroid/=area;
 
	return centroid;
 
}
 

	
 
/*Node &CellBase::getNode(list<Node *>::const_iterator i) const {
 

	
 
if (i==
 
	return m->getNode(i);
 
  }*/
 

	
 

	
 

	
 

	
 
void CellBase::SetIntegrals(void) const {
 
	
 
	// Set the initial values for the integrals over x^2,
 
	// xy, yy, x, and y
 
	
 
	// these values will be updated after each move of the CellBase wall
 
	
 
	intgrl_xx=0.; intgrl_xy=0.; intgrl_yy=0.;
 
	intgrl_x=0.; intgrl_y=0.;
 
	area=0.;
 
	list<Node *>::const_iterator nb;
 
	list<Node *>::const_iterator i=nodes.begin();
 
	
 
	for (;
 
		 i!=(nodes.end());
 
		 i++) {
 
		
 
		nb = i; nb++; if (nb==nodes.end()) nb=nodes.begin();
 
		
 
		area+=(*i)->x*(*nb)->y;
 
		area-=(*nb)->x*(*i)->y;
 
		intgrl_xx+= 
 
			((*i)->x*(*i)->x+
 
			 (*nb)->x*(*i)->x+
 
			 (*nb)->x*(*nb)->x ) *
 
			((*i)->x*(*nb)->y-
 
			 (*nb)->x*(*i)->y);
 
		intgrl_xy+= 
 
			((*nb)->x*(*i)->y-
 
			 (*i)->x*(*nb)->y)*
 
			((*i)->x*(2*(*i)->y+(*nb)->y)+
 
			 (*nb)->x*((*i)->y+2*(*nb)->y));
 
		intgrl_yy+=
 
			((*i)->x*(*nb)->y-
 
			 (*nb)->x*(*i)->y)*
 
			((*i)->y*(*i)->y+
 
			 (*nb)->y*(*i)->y+
 
			 (*nb)->y*(*nb)->y );
 
		intgrl_x+=
 
			((*nb)->x+(*i)->x)*
 
			((*i)->x*(*nb)->y-
 
			 (*nb)->x*(*i)->y);
 
		intgrl_y+=
 
			((*nb)->y+(*i)->y)*
 
			((*i)->x*(*nb)->y-
 
			 (*nb)->x*(*i)->y);
 
	}
 
	
 
	//area/=2.0;
 
	area = fabs(area)/2.0;
 
	
 
	/* intgrl_x/=6.;
 
	intgrl_y/=6.;
 
	
 
	intgrl_xx/=12.;
 
	intgrl_xy/=24.;
 
	intgrl_yy/=12.;*/
 
	
 
	
 
}
 

	
 
double CellBase::Length(Vector *long_axis, double *width)  const {
 
	
 
	// Calculate length and axes of CellBase
 
    
 
	// Calculate inertia tensor
 
	// see file inertiatensor.nb for explanation of this method
src/cellbase.h
Show inline comments
 
/*
 
 *
 
 *  $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;
 
	  if (division_axis) delete division_axis;
 
	  //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);
 
    void UnsetSource(void) {
 
      source = false;
 
    }
 
    
 
	inline bool Source(void) { return source; }
 
    enum boundary_type {None, Noflux, SourceSink, SAM};
 
    static const char * boundary_type_names[4];
 

	
 
    inline const char *BoundaryStr(void) { return boundary_type_names[boundary]; }
 
    
 
    ostream &print(ostream &os) const;
 

	
 
    inline double Chemical(int c) const { // returns the value of chemical c
 
    #ifdef _undefined_
 
      qDebug() << endl << "Entering cellbase::chemical(" << c << "), and nchem is: " << NChem() << "." << endl;
 
    #endif
 

	
 
      int nchem = NChem();
 

	
 
      #ifdef _undefined_
 
      if ((c<0) || (c>=nchem))
 
	MyWarning::warning("CellBase::Chemical says: index c is: %d, but nchem is: %d. Merely return zero", c, nchem);
 
      #endif
 

	
 
      return ((c<0) || (c>=nchem)) ? 0 : chem[c];
 
    }
 

	
 
    
 
    //void print_nblist(void) const;
 

	
 
    boundary_type SetBoundary(boundary_type bound) {
 
      if (bound!=None) {
 
	//area=0.;
 
	//length=0.;
 
      }
 
      return boundary=bound;
 
    }
 
  
 
    boundary_type ResetBoundary(void) {
 
      return boundary=None;
 
    }
 
    boundary_type Boundary(void) const {
 
      return boundary;
 
    }
 
    static int &NChem(void) {
 
      return static_data_members->nchem;
 
    }
 
  
 
    double CalcArea(void) const;
 
    double RecalcArea(void) {
 
      return area = CalcArea();
 
    }
 
    
 
    Vector Centroid(void) const;
 
  
 
    void SetIntegrals(void) const;
 
    
 
    double Length(Vector *long_axis = 0, double *width = 0) const;
 
    double CalcLength(Vector *long_axis = 0, double *width = 0) const;
 
    
 

	
 
    inline int Index(void) const {
 
      return index;
 
    }
 
  
 

	
 
    void SetTargetArea(double tar_ar) {
 
      target_area=tar_ar;
 
    }
 
    inline void SetTargetLength(double tar_l) {
 
      target_length=tar_l;
 
    }
 
    inline void SetLambdaLength(double lambda_length) {
 
      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 void DivideOverAxis(const Vector &v) {
 
		division_axis = new Vector(v);
 
		flag_for_divide = true;
 
	}
 
	
 
    //Vector Strain(void) const;
 

	
 
    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;
 
   	
 
	QString printednodelist(void);
 
 
 
   // void OnDivide(ParentInfo &parent_info, CellBase &daughter);
 

	
 
    
 
    inline bool DeadP(void) { return dead; }
 
    inline void MarkDead(void) { dead  = true; }
 
    
 
	static double &BaseArea(void) { 
 
		return static_data_members->base_area;
 
	}
 
  
 
    void CheckForDivision(void);
 

	
 
	
 
    // write flux from neighboring cells into "flux"
 
    void Flux(double *flux, double *D); 
 
    inline bool FixedP(void) { return fixed; }
 
    inline bool Fix(void) {  FixNodes(); return (fixed=true); }
 
    inline bool Unfix(void) { UnfixNodes(); return (fixed=false);}
 
    inline void setCellVec(Vector cv) { cellvec = cv; }
 
    
 
	bool AtBoundaryP(void) const;
 
    
 
    static inline int &NCells(void) {
 
      return static_data_members->ncells;
 
    }
 

	
 
 
 
   
 
    inline void Mark(void) {
 
      marked=true;
 
    }
 
    inline void Unmark(void) {
 
      marked=false;
 
    }
 
    inline bool Marked(void) const {
 
      return marked;
 
    }
 

	
 
	//! Returns the sum of chemical "chem" of this CellBase's neighbors
 
    double SumChemicalsOfNeighbors(int chem) {
 
      double sum=0.;
 
      for (list<Wall *>::const_iterator w=walls.begin();
 
	   w!=walls.end();
 
	   w++) {
 
	sum +=  (*w)->Length() * ( (*w)->c1!=this ? (*w)->c1->Chemical(chem) : (*w)->c2->Chemical(chem) );
 
      }
 
      return sum;
 
    }
 

	
 
    //! Generalization of the previous member function
 
    template<class P, class Op> P ReduceNeighbors(Op f) {
 
      P sum=0;
 
      for (list<Wall *>::const_iterator w=walls.begin();
 
	   w!=walls.end();
 
	   w++) {
 
	sum += (*w)->c1 != this ? f( *((*w)->c1) ) : f ( *((*w)->c2) ); 
 
      }
 
      return sum;
 
    }
 

	
 
    //! The same, but now for the walls
 
    template<class P, class Op> P ReduceWalls(Op f, P sum) {
 
      for (list<Wall *>::const_iterator w=walls.begin();
 
	   w!=walls.end();
 
	   w++) {
 
	sum += f( **w ); 
 
      }
 
      return sum;
 
    }
 
	
 
	
 
    
 
    
 
    //! The same, but now for the walls AND neighbors
 
    template<class P, class Op> P ReduceCellAndWalls(Op f) {
 
      P sum = 0;
 
      for (list<Wall *>::const_iterator w=walls.begin();
 
	   w!=walls.end();
 
	   w++) {
 
	sum += (*w)->c1 == this ? 
 
	  f( *((*w)->c1), *((*w)->c2), **w ) :  
 
	  f( *((*w)->c2), *((*w)->c1), **w );
 
      }
 
      return sum;
 
    }
 
    
 
	/* template<class Op> void LoopWalls(Op f) {
 
		for (list<Wall *>::const_iterator w=walls.begin();
 
			 w!=walls.end();
 
			 w++) {
 
			( **w)->f;
 
		}
 
	}*/
 
	
 
    //! Sum transporters at this CellBase's side of the walls
 
    double SumTransporters(int ch) {
 
      double sum=0.;
 
      for (list<Wall *>::const_iterator w=walls.begin();
 
	   w!=walls.end();
 
	   w++) {
 
	sum += (*w)->getTransporter(this, ch);
 
      
 
      }
 
      
 
      return sum;
 
    }
 

	
 
    inline int NumberOfDivisions(void) { return div_counter; }
 
    
 
    //! Sum transporters at this CellBase's side of the walls
 
    double SumLengthTransporters(int ch) {
 
      double sum=0.;
 
      for (list<Wall *>::const_iterator w=walls.begin();
 
	   w!=walls.end();
 
	   w++) {
 
	sum += (*w)->getTransporter(this, ch) * (*w)->Length();
 
      
 
      }
 
      
 
      return sum;
 
    }
 
    
 
	
 
    
 
    double SumLengthTransportersChemical(int trch, int ch) {
 
      double sum=0.;
 
      for (list<Wall *>::const_iterator w=walls.begin();
 
	   w!=walls.end();
 
	   w++) {
 
	sum += (*w)->getTransporter(this, trch) * ( (*w)->c1!=this ? (*w)->c1->Chemical(ch) : (*w)->c2->Chemical(ch) );
 
	
 
      }
 
      
 
      return sum;
 
    }
 
	inline int CellType(void) const { return cell_type; } 
 
	inline void SetCellType(int ct) { cell_type = ct; }
 

	
 
    
 
	static void SetNChem(int new_nchem) {
 
		if (NCells()) {
 
			MyWarning::error("CellBase::SetNChem says: not permitted, call SetNChem after deleting all cells.");
 
		} else {
 
			NChem() = new_nchem;
 
		}
 
	}
 
	
 
	inline double TargetLength() const { return target_length; } 
 

	
 
	static inline CellsStaticDatamembers *GetStaticDataMemberPointer(void) { return static_data_members; }
 
	
 
protected:
 
	// (define a list of Node* iterators)
 
	typedef list < list<Node *>::iterator > ItList;
 

	
 
    int index;
 

	
 
	inline void SetChemToNewchem(void) {
 
		for (int c=0;c<CellBase::NChem();c++) {
 
			chem[c]=new_chem[c];
 
		}
 
    }
 
    inline void SetNewChemToChem(void) {
 
		for (int c=0;c<CellBase::NChem();c++) {
 
			new_chem[c]=chem[c];
 
		}
 
    }
 
	inline double NewChem(int c) const { return new_chem[c]; }
 
	
 
 protected:
 
    list<Node *> nodes;
 
    void ConstructNeighborList(void);
 
	long wall_list_index (Wall *elem) const;
 

	
 
    // DATA MEMBERS
 
    
 
    // list of nodes, in clockwise order
 
  
 
    // 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;
 
	
 
	Vector *division_axis;
 
	
 
	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;
 
};
 

	
 
ostream &operator<<(ostream &os, const CellBase &v);
 

	
 
inline Vector PINdir(CellBase &here, CellBase &nb, Wall &w) {
 
	return w.getTransporter( &here, 1)  *  w.getInfluxVector(&here);
 
}
 

	
 

	
 
#endif
 

	
 

	
 

	
 

	
 

	
src/mesh.cpp
Show inline comments
 
@@ -1794,386 +1794,398 @@ protected:
 
    //cerr << "Calculated derivatives at " << x << "\n";    
 
  }
 
  
 
private:
 
  Mesh *m;
 
  int kmax,kount;
 
  double *xp,**yp,dxsav;
 
  bool monitor_window;
 
};
 
  
 

	
 

	
 
void Mesh::ReactDiffuse(double delta_t) {
 
  
 
  // Set Lengths of Walls
 
  for_each ( walls.begin(), walls.end(), 
 
	     mem_fun( &Wall::SetLength ) );
 
  
 
  static SolveMesh *solver = new SolveMesh(this);
 
  
 
  int nok, nbad, nvar;
 
  double *ystart = getValues(&nvar);
 
  
 
  solver->odeint(ystart, nvar, getTime(), getTime() + delta_t, 
 
		 par.ode_accuracy, par.dt, 1e-10, &nok, &nbad);
 
  
 
  setTime(getTime()+delta_t);
 
  setValues(getTime(),ystart);
 
  
 
}
 

	
 

	
 
Vector Mesh::FirstConcMoment(int chem) {
 
  
 
  Vector moment;
 
  for (vector<Cell *>::const_iterator c=cells.begin();
 
       c!=cells.end();
 
       c++) {
 
    
 
    moment += (*c)->Chemical(chem) * (*c)->Centroid();
 
    
 
  }
 

	
 
  return moment / (double)cells.size();
 
}
 

	
 
/*! This member function deletes all walls connected to two dead cells from the mesh.
 
  It should be called before the Cells are actually removed.
 
  If the cell is connect to one dead cell only, that reference is substituted for a reference 
 
  to the boundary polygon.
 
*/
 
void Mesh::DeleteLooseWalls(void) {
 

	
 
  list<Wall *>::iterator w=walls.begin();
 
  
 
  while (w!=walls.end()) {
 
    
 
    // if both cells of the wall are dead, remove the wall
 
    if ((*w)->C1()->DeadP() || (*w)->C2()->DeadP()) {
 
      if ((*w)->C1()->DeadP() && (*w)->C2()->DeadP()) {
 
	delete *w;
 
	w=walls.erase(w);
 
      } else {
 
	if ((*w)->C1()->DeadP())
 
	  (*w)->c1 = boundary_polygon;
 
	else
 
	  (*w)->c2 = boundary_polygon;
 
	w++;
 
      }
 
    } else {
 
      w++;
 
    }
 
    
 
  }
 
  
 
}
 

	
 
/*void Mesh::FitLeafToCanvas(double width, double height) {
 

	
 
  Vector bbll,bbur;
 
  BoundingBox(bbll,bbur);
 
  
 
  double scale_x = width/(bbur.x-bbll.x);
 
  double scale_y = height/(bbur.y-bbll.y);
 
  
 
  double factor = scale_x<scale_y ? scale_x:scale_y;
 
  
 
  Cell::SetMagnification(factor); // smallest of scale_x and scale_y
 
  
 
  double offset_x = (width/Cell::Magnification()-(bbur.x-bbll.x))/2.;  
 
  double offset_y = (height/Cell::Magnification()-(bbur.y-bbll.y))/2.;
 
  
 
  Cell::setOffset(offset_x, offset_y);
 
  
 
  }*/
 

	
 

	
 
void Mesh::CleanChemicals(const vector<double> &clean_chem, const vector<double> &clean_transporters) {
 
  
 
  if (clean_chem.size()!=(unsigned)Cell::NChem() || clean_transporters.size()!=(unsigned)Cell::NChem()) {
 
    throw "Run time error in Mesh::CleanChemicals: size of clean_chem and clean_transporters should be equal to Cell::NChem()";
 
  }
 
  for (vector<Cell *>::iterator c=cells.begin();
 
       c!=cells.end();
 
       c++) {
 
    
 
    for (int i=0;i<Cell::NChem();i++) {
 
      (*c)->SetChemical(i,clean_chem[i]);
 
    }
 
    (*c)->SetNewChemToChem();
 

	
 
  }
 

	
 
  // clean transporters
 
  for (list<Wall *>::iterator w=walls.begin();
 
       w!=walls.end();
 
       w++) {
 
    
 
    for (int i=0;i<Cell::NChem();i++) {
 
      (*w)->setTransporters1(i,clean_transporters[i]); (*w)->setNewTransporters1(i,clean_transporters[i]);
 
      (*w)->setTransporters2(i,clean_transporters[i]); (*w)->setNewTransporters2(i,clean_transporters[i]);
 
    }
 
  }
 
  
 
}
 

	
 
void Mesh::RandomizeChemicals(const vector<double> &max_chem, const vector<double> &max_transporters) {
 
  
 
  if (max_chem.size()!=(unsigned)Cell::NChem() || max_transporters.size()!=(unsigned)Cell::NChem()) {
 
    throw "Run time error in Mesh::CleanChemicals: size of max_chem and max_transporters should be equal to Cell::NChem()";
 
  }
 
  
 
  for (vector<Cell *>::iterator c=cells.begin();
 
       c!=cells.end();
 
       c++) {
 
    
 
    for (int i=0;i<Cell::NChem();i++) {
 
      (*c)->SetChemical(i,max_chem[i]*RANDOM());
 
    }
 
    (*c)->SetNewChemToChem();
 

	
 
  }
 

	
 
  // randomize transporters
 
  for (list<Wall *>::iterator w=walls.begin();
 
       w!=walls.end();
 
       w++) {
 
    
 
    for (int i=0;i<Cell::NChem();i++) {
 
      (*w)->setTransporters1(i,max_transporters[i] * RANDOM()); (*w)->setNewTransporters1(i, (*w)->Transporters1(i) );
 
      (*w)->setTransporters2(i,max_transporters[i] * RANDOM()); (*w)->setNewTransporters2(i, (*w)->Transporters1(i) );
 
    }
 
  }
 
  
 
}
 

	
 
//!\brief Calculates a vector with derivatives of all variables, which
 
// we can pass to an ODESolver. 
 
void Mesh::Derivatives(double *derivs) {
 
  
 
  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;
 
  
 
  //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
 
	  plugin->CelltoCellTransport(*w, &(derivs[(*w)->c1->Index() * nchems]),
 
								  &(derivs[(*w)->c2->Index() * nchems]));
 
	  double *dchem_c1 = &(derivs[(*w)->c1->Index() * nchems]);
 
	  double *dchem_c2 = &(derivs[(*w)->c2->Index() * nchems]);
 
	  //plugin->CelltoCellTransport(*w, &(derivs[(*w)->c1->Index() * nchems]),
 
							//	  &(derivs[(*w)->c2->Index() * nchems]));
 
	  // quick fix: dummy values to prevent end user from writing into outer space and causing a crash :-)
 
	  // start here if you want to implement chemical input/output into environment over boundaries
 
	  double dummy1, dummy2;
 
	  if ((*w)->c1->Index()<0) { // tests if c1 is the boundary pol
 
		  dchem_c1 = &dummy1;
 
	  }
 
	  if ((*w)->c2->Index()<0) {
 
		  dchem_c2 = &dummy2;
 
	  }
 
	  plugin->CelltoCellTransport(*w, dchem_c1, dchem_c2); 
 
	  
 
	  //(*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 c=cells.begin();
 
       c!=cells.end();
 
       c++) {
 
    for (int ch=0;ch<nchems;ch++) {
 
      (*c)->SetChemical(ch, y[i+ch]);
 
    }
 
    if ( !(emit_count%stride)) {
 
      (*c)->EmitValues(x);
 
    }
 
    i+=nchems;
 
  }
 
  
 
  for (list<Wall *>::iterator w=walls.begin();
 
       w!=walls.end();
 
       w++) {
 
    
 
    for (int ch=0;ch<nchems;ch++) {
 
      (*w)->setTransporters1(ch,y[i+ch]);
 
    }
 
    i+=nchems;
 
    
 
    for (int ch=0;ch<nchems;ch++) {
 
      (*w)->setTransporters2(ch,y[i+ch]);
 
    }
 
    i+=nchems;
 

	
 
  }
 

	
 
  emit_count++;
 
}
 

	
 
double *Mesh::getValues(int *neqs) {
 
  
 
  int nwalls = walls.size();
 
  int ncells = cells.size();
 
  int nchems = Cell::NChem();
 
  
 
  // two eqs per chemical for each wall, 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.
 
  (*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) ] ]
 

	
 
  static double *values = 0;
 
  if (values!=0) { delete[] values; }
 
  
 
  values = new double[*neqs];
 
  
 
  int i=0;
 
  for (vector<Cell *>::iterator c=cells.begin();
 
       c!=cells.end();
 
       c++) {
 
    for (int ch=0;ch<nchems;ch++) {
 
      values[i+ch]=(*c)->Chemical(ch);
 
    }
 
    i+=nchems;
 
  }
 
  
 
  for (list<Wall *>::iterator w=walls.begin();
 
       w!=walls.end();
 
       w++) {
 
    
 
    for (int ch=0;ch<nchems;ch++) {
 
      values[i+ch]=(*w)->Transporters1(ch);
 
    }
 
    i+=nchems;
 
    
 
    for (int ch=0;ch<nchems;ch++) {
 
      values[i+ch]=(*w)->Transporters2(ch);
 
    }
 
    i+=nchems;
 

	
 
  }
 

	
 
  return values;
 
}
 

	
 
void Mesh::DrawNodes(QGraphicsScene *c) const {
 
  
 
  for (vector<Node *>::const_iterator n=nodes.begin();
 
       n!=nodes.end();
 
       n++) {
 
    
 
    Node *i=*n;
 
    
 
    NodeItem *item = new NodeItem ( &(*i), c );
 
    item->setColor();
 
    
 
    item->setZValue(5);
 
    item->show();
 
    item ->setPos(((Cell::offset[0]+i->x)*Cell::factor),
 
		  ((Cell::offset[1]+i->y)*Cell::factor) );
 
  }
 

	
 
}
 

	
 
/*! Returns the sum of protein "ch" of a cycling protein in cells and walls */
 
double Mesh::CalcProtCellsWalls(int ch) const {
 

	
 

	
 
  double sum_prot=0.;
 

	
 
  // At membranes
 
  for (list<Wall *>::const_iterator w=walls.begin();
 
       w!=walls.end();
 
       w++) {
 
    sum_prot += (*w)->Transporters1(ch);
 
    sum_prot += (*w)->Transporters2(ch);
 
  }
 

	
 
  // At cells
 
  for (vector<Cell *>::const_iterator c=cells.begin();
 
       c!=cells.end();
 
       c++) {
 
    
 
    sum_prot += (*c)->Chemical(ch);
 
  }
 
  
 
  return sum_prot;
 

	
 
}
 

	
 
void Mesh::SettoInitVals(void) {
 
  
 
  vector<double> clean_chem(Cell::NChem());
 
  vector<double> clean_transporters(Cell::NChem());
 
  
 
  for (int i=0;i<Cell::NChem();i++) {
 
    clean_transporters[i]=0.;
 
    clean_chem[i]=par.initval[i];
 
  }
 

	
 
  // Amount of PIN1
 
  //clean_chem[1] = 0.;
 
  
 
  CleanChemicals(clean_chem, clean_transporters);
 

	
 
}
 

	
 
string Mesh::getTimeHours(void) const {
 
	int hours = (int)(time / 3600);
 
	int mins = (int)((time - hours * 3600)/60);
 
	int secs = (int)((time - hours * 3600 - mins * 60));
 
	ostringstream tstr;
 
	tstr << hours << " h " << mins << " m " << secs << " s";
 
	return tstr.str();
 
}
 

	
 
QVector<qreal> Mesh::VertexAngles(void) {
 
	QVector<qreal> angles;
 
	for (vector<Node *>::const_iterator n=nodes.begin();
 
		 n!=nodes.end();
 
		 n++) {
 
		if ((*n)->Value()>2 && !(*n)->BoundaryP() ) {
 
			angles+=(*n)->NeighbourAngles();
 
		}
 
	}
 
	return angles;
 
}
 

	
 
QVector< QPair<qreal,int> > Mesh::VertexAnglesValues(void) {
 

	
 
	QVector< QPair<qreal,int> > anglesvalues;
 
	for (vector<Node *>::const_iterator n=nodes.begin();
 
		 n!=nodes.end();
 
		 n++) {
src/mesh.h
Show inline comments
 
@@ -53,385 +53,391 @@ public:
 
			queue<T,C>::c.push_back(x);
 
		}
 
	}
 
	void clear(void) {
 
		queue<T,C>::c.clear();
 
	}
 
};
 

	
 
//template<class P> P& deref_ptr<P>( P *obj) { return *obj; }
 
template<class P> P& deref_ptr ( P *obj) { return *obj; }
 

	
 

	
 
class Mesh {
 
	
 
	friend class Cell;
 
	friend class Node;
 
	friend class FigureEditor;
 
	
 
public: 
 
	Mesh(void) {
 
		// Make sure the reserved value is large enough if a cell is added
 
		// in "Divide" when the capacity is insufficient, "cells" might be
 
		// relocated including the current Cell (i.e. the value of *this)
 
		// calling "Mesh::IncreaseCapacityIfNecessary" (from another
 
		// object than Cell, e.g. Mesh) before entering Divide will solve
 
		// this issue (solved now).
 
		cells.reserve(2);
 
		nodes.reserve(500);
 
		
 
		//boundary_polygon = new BoundaryPolygon();
 
		
 
		time = 0.;
 
		plugin = 0;
 
	};
 
	~Mesh(void) {
 
		delete boundary_polygon;
 
		/* if (plugin)
 
			delete plugin;*/
 
	};
 
	
 
	void Clean(void);
 
	//void Plane(int xwidth, int ywidth, int nx, int ny, bool randomP=false);
 
	Cell &EllipticCell(double xc, double yc, double ra, double rb, int nnodes=10, double rotation=0);
 
	Cell &CircularCell(double xc, double yc, double r, int nnodes=10) {
 
		return EllipticCell(xc, yc, r, r, nnodes, 0);
 
	}
 
	Cell &LeafPrimordium(int n, double pet_length);
 
	Cell &LeafPrimordium2(int n);
 
	Cell *RectangularCell(const Vector ll, const Vector ur, double rotation = 0);
 
	void CellFiles(const Vector ll, const Vector ur);
 
	
 
	/*  void GMVoutput(ostream &os, 
 
	 const char *codename=0, const char *codever=0,
 
	 const char *comments=0);*/
 
	//void RandPoints(int npoints);
 
	
 
	inline Cell &getCell(int i) {
 
		if ((unsigned)i<cells.size())
 
			return *cells[i];
 
		else {
 
			cerr << i << endl;
 
			cerr << "size is " << cells.size() << endl;
 
			abort();
 
			//	throw("Index out of range in Mesh::getCell");
 
		}
 
	}
 
	
 
	inline Node &getNode(int i) {
 
		//if (i >= nodes.size() || i < 0) {
 
		//  cerr << "Mesh::getNode: Warning. Index " << i << " out of range.\n";
 
		// }
 
		return *nodes[i];    
 
	}
 
	
 
	//double Diffusion(void);
 
	inline int size(void) {
 
		return cells.size();
 
	}
 
	inline int nnodes(void) {
 
		return nodes.size();
 
	}
 
	//void SortNBLists(void);
 
	
 
	/*template<class Op> void LoopCells(Op f) {
 
	 for (vector<Cell>::iterator i=cells.begin();
 
	 i!=cells.end();
 
	 i++) {
 
	 f(*i); 
 
	 }
 
	 }*/
 
	
 
	/*! \brief Calls function f for all Cells f.
 
	 
 
	 Using this template requires some fiddling with function adaptors bind2nd and mem_fun_ref.
 
	 
 
	 Example usage for calling a member function on each cell:
 
	 
 
	 mesh.LoopCells( bind2nd (mem_fun_ref( &Cell::DrawDiffEdges), &canvas ) );
 
	 
 
	 This calls Cell's member function DrawDiffEdges, taking one
 
	 argument canvas, on all Cell objects in the current Mesh.
 
	 
 
	 */
 
	template<class Op> void LoopCells(Op f) {
 
		for (vector <Cell *>::iterator i=cells.begin();
 
			 i!=cells.end();
 
			 i++) {
 
			f(**i);
 
		}
 
		//for_each(cells.begin(),cells.end(),f);
 
	}
 
	
 
	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);
 
			
 
		}
 
		//for_each(cells.begin(),cells.end(),f);
 
	}
 
	
 
	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++) {
 
			f(*shuffled_nodes[*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++) {
 
			f(*shuffled_cells[*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;
 
			}
 
		}
 
	}
 
/*	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();
 
		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;
 
		
 
		return neqs;
 
	}
 
	void IncreaseCellCapacityIfNecessary(void) {
 
		
 
		return;
 
		// cerr << "Entering Mesh::IncreaseCellCapacityIfNecessary \n";
 
		// make sure we always have enough space 
 
		// to have each cell divide at least once
 
		//
 
		// Note that we must do this, because Cell::Divide pushes a new Cell
 
		// onto Mesh::cells. As a result, Mesh::cells might be relocated 
 
		// if we are _within_ a Cell object: i.e. pointer "this" will be changed!!
 
		// 
 
		// An alternative solution could be to make "Mesh::cells" a list,
 
		// but this won't work because we need random access for 
 
		// the Monte Carlo algorithm.
 
		
 
		if (2*cells.size()>cells.capacity()) {
 
			cerr << "Increasing capacity to "  << 2*cells.capacity() << endl;
 
			cerr << "Current capacity is " << cells.capacity() << endl;
 
			cells.reserve(cells.capacity()*2);
 
		}
 
	}
 
	
 
	void ReserveMoreCells(int n) {
 
		if (nodes.size()+n>nodes.capacity()) {
 
			nodes.reserve(size()+n);
 
		}
 
	}
 
	double Area(void);
 
	double MeanArea(void) {
 
		double sum=0.;
 
		for (vector<Cell *>::const_iterator i=cells.begin();
 
			 i!=cells.end();
 
			 i++) {
 
			sum+=(*i)->Area();
 
		}
 
		return sum/(double)NCells();
 
	}
 
	
 
	void SetBaseArea(void);
 
	int NCells(void) const {
 
		return cells.size();
 
	}
 
	inline int NNodes(void) const {
 
		return nodes.size();
 
	}
 
	void PrintQueue(ostream &os) {
 
		while (!node_insertion_queue.empty()) {
 
			os << node_insertion_queue.front() << endl;
 
			node_insertion_queue.pop();
 
		}
 
		//copy (node_insertion_queue.begin(),node_insertion_queue.end(),ostream_iterator<Edge>(cerr, " "));
 
	}
 
	
 
	void InsertNodes(void) {
 
		// insert the nodes in the insertion queue
 
		while (!node_insertion_queue.empty()) {
 
			
 
			//cerr << node_insertion_queue.front() << endl;
 
			InsertNode(node_insertion_queue.front());
 
			node_insertion_queue.pop();
 
		}
 
		
 
	}
 
	
 
	void Clear(); 
 
	
 
	//  template<class ReactFunction> ReactDiffuse(const double D, ReactFunction& react) {
 
	void ReactDiffuse( double delta_t = 1 );
 
	//void Diffuse(const double D);
 
	double SumChemical(int ch);
 
	void SetChemical(int ch, double value) {
 
		for (vector<Cell *>::iterator c=cells.begin();
 
			 c!=cells.end();
 
			 c++) {
 
			(*c)->chem[ch]=value;
 
		}
 
	}
 
	
 
	// used for interacing with ODE-solvers (e.g. NRCRungeKutta)
 
	void setValues(double x, double *y);
 
	double *getValues(int *neqs);
 
	void Derivatives(double *derivs);
 
#ifdef QTGRAPHICS
 
	inline void DrawBoundary(QGraphicsScene *c) {
 
		boundary_polygon->Draw(c);
 
	}
 
	void DrawNodes(QGraphicsScene *c) const;
 
	
 
#endif
 
	double max_chem;
 
	
 
	void XMLSave(const char *docname, xmlNode *settings=0) const;
 
	void XMLRead(const char *docname, xmlNode **settings=0, bool geometry = true, bool pars = true, bool simtime = true);
 
	void XMLReadPars(const xmlNode * root_node);
 
	void XMLReadGeometry(const xmlNode *root_node);
 
	void XMLReadSimtime(const xmlNode *root_node);
 
	void XMLReadNodes(xmlNode *cur);
 
	void XMLReadCells(xmlNode *cur);
 
	void XMLParseTree(const xmlNode * root_node);
 
	void XMLReadWalls(xmlNode *cur, vector<Wall *> *tmp_cells);
 
	void XMLReadWallsToCells(xmlNode *root, vector<Wall *> *tmp_walls);
 
	void XMLReadNodeSets(xmlNode *root);
 
	void XMLReadNodeSetsToNodes(xmlNode *root);
 
	void PerturbChem(int chemnum, double range);
 
	void CleanUpCellNodeLists(void);
 
	void CleanUpWalls(void);
 
	void CutAwayBelowLine( Vector startpoint, Vector endpoint );
 
	void CutAwaySAM(void);
 
	void RepairBoundaryPolygon(void);
 
	void Rotate(double angle, Vector center);
 
	void PrintWallList( void );
 
	void TestIllegalWalls(void);
 
	Vector FirstConcMoment(int chem);
 
	inline Vector Centroid(void) {
 
		return boundary_polygon->Centroid();
 
	}
 
	
 
	inline Vector Offset(void) {
 
		return boundary_polygon->Offset();
 
	}
 
	
 
	inline double Factor(void) {
 
		return boundary_polygon->Factor();
 
	}
 
	
 
	void DeleteLooseWalls(void);
 
	void FitLeafToCanvas(double width, double height);
 
	void AddNodeSet(NodeSet *node_set) {
 
		node_sets.push_back(node_set);
 
	}
 
	
 
	void CleanChemicals(const vector<double> &clean_chem, const vector<double> &clean_transporters);
 
	void RandomizeChemicals(const vector<double> &max_chem, const vector<double> &max_transporters);
 
	inline double getTime(void) const { return time; }
 
	string getTimeHours(void) const; 
 
	inline void setTime(double t) { time = t; }
 
	double CalcProtCellsWalls(int ch) const;  
 
	void SettoInitVals(void);
 
	QVector<qreal> VertexAngles(void);
 
	QVector< QPair<qreal,int> > VertexAnglesValues(void);
 
	void SetSimPlugin(SimPluginInterface *new_plugin) {
 
		/* if (plugin) 
 
			delete plugin;*/
 
		plugin=new_plugin;
 
	}
 
	QString ModelID(void) { return plugin?plugin->ModelID():QString("undefined"); }
 
	void StandardInit(void);	
 

	
 
	Node* findNextBoundaryNode(Node*);
 

	
 
private:
 
	
 
	// Data members
 
	vector<Cell *> cells;
 
	vector<Node *> nodes;
 
	list<Wall *> walls; // we need to erase elements from this container frequently, hence a list.
 
public:
 
	vector<NodeSet *> node_sets;
 
private:
0 comments (0 inline, 0 general)