EI/VirtualLeaf Changeset - 413a4067bfff · Centrum Wiskunde & Informatica (CWI)

Changeset - 413a4067bfff

Parent rev.

Child rev.

[Not reviewed]

default

0 17 19

Roeland Merks - 15 years ago 2010-06-03 18:40:26
roeland.merks@cwi.nl

Changed interface for model plugins to make usage of pointers and references consistent.
Added codes for tutorial belonging to VLeaf manuscript.
Something changed in modelcatalogue. I'm not sure what...

user: Roeland Merks <roeland.merks@cwi.nl>
branch 'default'
added src/TutorialCode/Tutorial0/mymodel.cpp
added src/TutorialCode/Tutorial0/mymodel.h
added src/TutorialCode/Tutorial0/mymodel.pro
added src/TutorialCode/Tutorial1/mymodel.cpp
added src/TutorialCode/Tutorial1/mymodel.h
added src/TutorialCode/Tutorial1/mymodel.pro
added src/TutorialCode/Tutorial2/mymodel.cpp
added src/TutorialCode/Tutorial2/mymodel.h
added src/TutorialCode/Tutorial2/mymodel.pro
added src/TutorialCode/Tutorial3/mymodel.cpp
added src/TutorialCode/Tutorial3/mymodel.h
added src/TutorialCode/Tutorial3/mymodel.pro
added src/TutorialCode/Tutorial4/mymodel.cpp
added src/TutorialCode/Tutorial4/mymodel.h
added src/TutorialCode/Tutorial4/mymodel.pro
added src/TutorialCode/Tutorial5/mymodel.cpp
added src/TutorialCode/Tutorial5/mymodel.h
added src/TutorialCode/Tutorial5/mymodel.pro
added src/build_models/translate_plugin.pl
changed src/VirtualLeafpar.tmpl
changed src/build_models/auxingrowthplugin.cpp
changed src/build_models/auxingrowthplugin.h
changed src/build_models/meinhardtplugin.cpp
changed src/build_models/meinhardtplugin.h
changed src/build_models/plugin_auxingrowth.pro
changed src/build_models/plugin_leaf.pro
changed src/build_models/plugin_meinhardt.pro
changed src/build_models/plugin_test.pro
changed src/cell.cpp
changed src/cellbase.cpp
changed src/cellbase.h
changed src/libplugin.pro
changed src/mesh.cpp
changed src/mesh.h
changed src/modelcatalogue.cpp
changed src/simplugin.h

36 files changed with 1397 insertions and 75 deletions:

src/TutorialCode/Tutorial0/mymodel.cpp

src/TutorialCode/Tutorial0/mymodel.h

src/TutorialCode/Tutorial0/mymodel.pro

src/TutorialCode/Tutorial1/mymodel.cpp

src/TutorialCode/Tutorial1/mymodel.h

src/TutorialCode/Tutorial1/mymodel.pro

src/TutorialCode/Tutorial2/mymodel.cpp

src/TutorialCode/Tutorial2/mymodel.h

src/TutorialCode/Tutorial2/mymodel.pro

src/TutorialCode/Tutorial3/mymodel.cpp

src/TutorialCode/Tutorial3/mymodel.h

src/TutorialCode/Tutorial3/mymodel.pro

src/TutorialCode/Tutorial4/mymodel.cpp

src/TutorialCode/Tutorial4/mymodel.h

src/TutorialCode/Tutorial4/mymodel.pro

src/TutorialCode/Tutorial5/mymodel.cpp

120

src/TutorialCode/Tutorial5/mymodel.h

src/TutorialCode/Tutorial5/mymodel.pro

src/VirtualLeafpar.tmpl

src/build_models/auxingrowthplugin.cpp

src/build_models/auxingrowthplugin.h

src/build_models/meinhardtplugin.cpp

src/build_models/meinhardtplugin.h

src/build_models/plugin_auxingrowth.pro

src/build_models/plugin_leaf.pro

src/build_models/plugin_meinhardt.pro

src/build_models/plugin_test.pro

src/build_models/translate_plugin.pl

src/cell.cpp

src/cellbase.cpp

src/cellbase.h

src/libplugin.pro

src/mesh.cpp

src/mesh.h

src/modelcatalogue.cpp

src/simplugin.h

0 comments (0 inline, 0 general)

src/TutorialCode/Tutorial0/mymodel.cpp

➞

Show inline comments

@@ new file 100644 @@
 /*
+ *
  *  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 <QObject>
 #include <QtGui>
 #include "simplugin.h"
 #include "parameter.h"
 #include "wallbase.h"
 #include "cellbase.h"
 #include "mymodel.h"
 static const std::string _module_id("$Id$");
 QString MyModel::ModelID(void) {
   // specify the name of your model here
   return QString( "My first model in VirtualLeaf" );
+}
 // return the number of chemicals your model uses
 int MyModel::NChem(void) { return 0; }
 // To be executed after cell division
 void MyModel::OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2) {
   // rules to be executed after cell division go here
   // (e.g., cell differentiation rules)
+}
 void MyModel::SetCellColor(CellBase *c, QColor *color) {
   // add cell coloring rules here
+}
 void MyModel::CellHouseKeeping(CellBase *c) {
   // add cell behavioral rules here
+}
 void MyModel::CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2) {
   // add biochemical transport rules here
+}
 void MyModel::WallDynamics(Wall *w, double *dw1, double *dw2) {
   // add biochemical networks for reactions occuring at walls here
+}
 void MyModel::CellDynamics(CellBase *c, double *dchem) {
   // add biochemical networks for intracellular reactions here
+}
 Q_EXPORT_PLUGIN2(mymodel, MyModel)

src/TutorialCode/Tutorial0/mymodel.h

➞

Show inline comments

@@ new file 100644 @@
 /*
  *  $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.
+ *
  */
 #include <QObject>
 #include <QtGui>
 #include <QString>
 #include "simplugin.h"
 class MyModel : public QObject, SimPluginInterface {
 	Q_OBJECT
 	Q_INTERFACES(SimPluginInterface);
 public:
 	virtual QString ModelID(void);
 	// Executed after the cellular mechanics steps have equillibrized
 	virtual void CellHouseKeeping (CellBase *c);
 	// Differential equations describing transport of chemicals from cell to cell
 	virtual void CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2);
 	// Differential equations describing chemical reactions taking place at or near the cell walls
 	// (e.g. PIN accumulation)
 	virtual void WallDynamics(Wall *w, double *dw1, double *dw2);
 	// Differential equations describing chemical reactions inside the cells
 	virtual void CellDynamics(CellBase *c, double *dchem);
 	// to be executed after a cell division
 	virtual void OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2);
 	// to be executed for coloring a cell
 	virtual void SetCellColor(CellBase *c, QColor *color);
 	// return number of chemicals
 	virtual int NChem(void);
 };

src/TutorialCode/Tutorial0/mymodel.pro

➞

Show inline comments

@@ new file 100644 @@
+#
 # $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.
+#
 TARGET = mymodel
 VLEAFHOME = ../../..
 CONFIG += release
 CONFIG -= debug
 CONFIG += plugin
 BINDIR = $${VLEAFHOME}/bin
 LIBDIR = $${VLEAFHOME}/lib
 INCDIR = $${VLEAFHOME}/src
 DEFINES = QTGRAPHICS # VLEAFPLUGIN
 DESTDIR = $${BINDIR}/models
 HEADERS = $${TARGET}.h
 INCLUDEPATH += $${INCDIR}
 QMAKE_CXXFLAGS += -fexceptions #-I$${INCDIR}
 QMAKE_CXXFLAGS_DEBUG += -g3
 QMAKE_CXXFLAGS_DEBUG += -DQDEBUG
 QT += qt3support
 SOURCES = $${TARGET}.cpp
 TEMPLATE = lib
 unix {
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2
  QMAKE_LFLAGS += -fPIC
+}
 win32 {
  LIBXML2DIR = C:\libxml2
  LIBICONVDIR = C:\libiconv
  LIBZDIR = C:\libz
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -DLIBXML_STATIC
  QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include
+}
 # finish

src/TutorialCode/Tutorial1/mymodel.cpp

➞

Show inline comments

@@ new file 100644 @@
 /*
+ *
  *  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 <QObject>
 #include <QtGui>
 #include "simplugin.h"
 #include "parameter.h"
 #include "wallbase.h"
 #include "cellbase.h"
 #include "mymodel.h"
 static const std::string _module_id("$Id$");
 QString MyModel::ModelID(void) {
   // specify the name of your model here
   return QString( "Cell growth" );
+}
 // return the number of chemicals your model uses
 int MyModel::NChem(void) { return 0; }
 // To be executed after cell division
 void MyModel::OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2) {
   // rules to be executed after cell division go here
   // (e.g., cell differentiation rules)
+}
 void MyModel::SetCellColor(CellBase *c, QColor *color) {
   // add cell coloring rules here
+}
 void MyModel::CellHouseKeeping(CellBase *c) {
   // add cell behavioral rules here
 	c->EnlargeTargetArea(par->cell_expansion_rate);
+}
 void MyModel::CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2) {
   // add biochemical transport rules here
+}
 void MyModel::WallDynamics(Wall *w, double *dw1, double *dw2) {
   // add biochemical networks for reactions occuring at walls here
+}
 void MyModel::CellDynamics(CellBase *c, double *dchem) {
   // add biochemical networks for intracellular reactions here
+}
 Q_EXPORT_PLUGIN2(mymodel, MyModel)

src/TutorialCode/Tutorial1/mymodel.h

➞

Show inline comments

@@ new file 100644 @@
 /*
  *  $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.
+ *
  */
 #include <QObject>
 #include <QtGui>
 #include <QString>
 #include "simplugin.h"
 class MyModel : public QObject, SimPluginInterface {
 	Q_OBJECT
 	Q_INTERFACES(SimPluginInterface);
 public:
 	virtual QString ModelID(void);
 	// Executed after the cellular mechanics steps have equillibrized
 	virtual void CellHouseKeeping (CellBase *c);
 	// Differential equations describing transport of chemicals from cell to cell
 	virtual void CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2);
 	// Differential equations describing chemical reactions taking place at or near the cell walls
 	// (e.g. PIN accumulation)
 	virtual void WallDynamics(Wall *w, double *dw1, double *dw2);
 	// Differential equations describing chemical reactions inside the cells
 	virtual void CellDynamics(CellBase *c, double *dchem);
 	// to be executed after a cell division
 	virtual void OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2);
 	// to be executed for coloring a cell
 	virtual void SetCellColor(CellBase *c, QColor *color);
 	// return number of chemicals
 	virtual int NChem(void);
 };

src/TutorialCode/Tutorial1/mymodel.pro

➞

Show inline comments

@@ new file 100644 @@
+#
 # $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.
+#
 TARGET = mymodel
 VLEAFHOME = ../../..
 CONFIG += release
 CONFIG -= debug
 CONFIG += plugin
 BINDIR = $${VLEAFHOME}/bin
 LIBDIR = $${VLEAFHOME}/lib
 INCDIR = $${VLEAFHOME}/src
 DEFINES = QTGRAPHICS # VLEAFPLUGIN
 DESTDIR = $${BINDIR}/models
 HEADERS = $${TARGET}.h
 INCLUDEPATH += $${INCDIR}
 QMAKE_CXXFLAGS += -fexceptions #-I$${INCDIR}
 QMAKE_CXXFLAGS_DEBUG += -g3
 QMAKE_CXXFLAGS_DEBUG += -DQDEBUG
 QT += qt3support
 SOURCES = $${TARGET}.cpp
 TEMPLATE = lib
 unix {
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2
  QMAKE_LFLAGS += -fPIC
+}
 win32 {
  LIBXML2DIR = C:\libxml2
  LIBICONVDIR = C:\libiconv
  LIBZDIR = C:\libz
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -DLIBXML_STATIC
  QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include
+}
 # finish

src/TutorialCode/Tutorial2/mymodel.cpp

➞

Show inline comments

@@ new file 100644 @@
 /*
+ *
  *  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 <QObject>
 #include <QtGui>
 #include "simplugin.h"
 #include "parameter.h"
 #include "wallbase.h"
 #include "cellbase.h"
 #include "mymodel.h"
 static const std::string _module_id("$Id$");
 QString MyModel::ModelID(void) {
   // specify the name of your model here
   return QString( "Cell growth and cell division" );
+}
 // return the number of chemicals your model uses
 int MyModel::NChem(void) { return 0; }
 // To be executed after cell division
 void MyModel::OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2) {
   // rules to be executed after cell division go here
   // (e.g., cell differentiation rules)
+}
 void MyModel::SetCellColor(CellBase *c, QColor *color) {
   // add cell coloring rules here
+}
 void MyModel::CellHouseKeeping(CellBase *c) {
   // add cell behavioral rules here
 	c->EnlargeTargetArea(par->cell_expansion_rate);
 	if (c->Area() > par->rel_cell_div_threshold * c->BaseArea()) {
 		c->Divide();
+	}
+}
 void MyModel::CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2) {
   // add biochemical transport rules here
+}
 void MyModel::WallDynamics(Wall *w, double *dw1, double *dw2) {
   // add biochemical networks for reactions occuring at walls here
+}
 void MyModel::CellDynamics(CellBase *c, double *dchem) {
   // add biochemical networks for intracellular reactions here
+}
 Q_EXPORT_PLUGIN2(mymodel, MyModel)

src/TutorialCode/Tutorial2/mymodel.h

➞

Show inline comments

@@ new file 100644 @@
 /*
  *  $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.
+ *
  */
 #include <QObject>
 #include <QtGui>
 #include <QString>
 #include "simplugin.h"
 class MyModel : public QObject, SimPluginInterface {
 	Q_OBJECT
 	Q_INTERFACES(SimPluginInterface);
 public:
 	virtual QString ModelID(void);
 	// Executed after the cellular mechanics steps have equillibrized
 	virtual void CellHouseKeeping (CellBase *c);
 	// Differential equations describing transport of chemicals from cell to cell
 	virtual void CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2);
 	// Differential equations describing chemical reactions taking place at or near the cell walls
 	// (e.g. PIN accumulation)
 	virtual void WallDynamics(Wall *w, double *dw1, double *dw2);
 	// Differential equations describing chemical reactions inside the cells
 	virtual void CellDynamics(CellBase *c, double *dchem);
 	// to be executed after a cell division
 	virtual void OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2);
 	// to be executed for coloring a cell
 	virtual void SetCellColor(CellBase *c, QColor *color);
 	// return number of chemicals
 	virtual int NChem(void);
 };

src/TutorialCode/Tutorial2/mymodel.pro

➞

Show inline comments

@@ new file 100644 @@
+#
 # $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.
+#
 TARGET = mymodel
 VLEAFHOME = ../../..
 CONFIG += release
 CONFIG -= debug
 CONFIG += plugin
 BINDIR = $${VLEAFHOME}/bin
 LIBDIR = $${VLEAFHOME}/lib
 INCDIR = $${VLEAFHOME}/src
 DEFINES = QTGRAPHICS # VLEAFPLUGIN
 DESTDIR = $${BINDIR}/models
 HEADERS = $${TARGET}.h
 INCLUDEPATH += $${INCDIR}
 QMAKE_CXXFLAGS += -fexceptions #-I$${INCDIR}
 QMAKE_CXXFLAGS_DEBUG += -g3
 QMAKE_CXXFLAGS_DEBUG += -DQDEBUG
 QT += qt3support
 SOURCES = $${TARGET}.cpp
 TEMPLATE = lib
 unix {
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2
  QMAKE_LFLAGS += -fPIC
+}
 win32 {
  LIBXML2DIR = C:\libxml2
  LIBICONVDIR = C:\libiconv
  LIBZDIR = C:\libz
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -DLIBXML_STATIC
  QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include
+}
 # finish

src/TutorialCode/Tutorial3/mymodel.cpp

➞

Show inline comments

@@ new file 100644 @@
 /*
+ *
  *  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 <QObject>
 #include <QtGui>
 #include "simplugin.h"
 #include "parameter.h"
 #include "wallbase.h"
 #include "cellbase.h"
 #include "mymodel.h"
 static const std::string _module_id("$Id$");
 QString MyModel::ModelID(void) {
   // specify the name of your model here
   return QString( "Cell growth and cell division" );
+}
 // return the number of chemicals your model uses
 int MyModel::NChem(void) { return 0; }
 // To be executed after cell division
 void MyModel::OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2) {
   // rules to be executed after cell division go here
   // (e.g., cell differentiation rules)
+}
 void MyModel::SetCellColor(CellBase *c, QColor *color) {
   // add cell coloring rules here
+}
 void MyModel::CellHouseKeeping(CellBase *c) {
   // add cell behavioral rules here
 	c->EnlargeTargetArea(par->cell_expansion_rate);
 	if (c->Area() > par->rel_cell_div_threshold * c->BaseArea()) {
 		c->DivideOverAxis(Vector(1,0,0));
+	}
+}
 void MyModel::CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2) {
   // add biochemical transport rules here
+}
 void MyModel::WallDynamics(Wall *w, double *dw1, double *dw2) {
   // add biochemical networks for reactions occuring at walls here
+}
 void MyModel::CellDynamics(CellBase *c, double *dchem) {
   // add biochemical networks for intracellular reactions here
+}
 Q_EXPORT_PLUGIN2(mymodel, MyModel)

src/TutorialCode/Tutorial3/mymodel.h

➞

Show inline comments

@@ new file 100644 @@
 /*
  *  $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.
+ *
  */
 #include <QObject>
 #include <QtGui>
 #include <QString>
 #include "simplugin.h"
 class MyModel : public QObject, SimPluginInterface {
 	Q_OBJECT
 	Q_INTERFACES(SimPluginInterface);
 public:
 	virtual QString ModelID(void);
 	// Executed after the cellular mechanics steps have equillibrized
 	virtual void CellHouseKeeping (CellBase *c);
 	// Differential equations describing transport of chemicals from cell to cell
 	virtual void CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2);
 	// Differential equations describing chemical reactions taking place at or near the cell walls
 	// (e.g. PIN accumulation)
 	virtual void WallDynamics(Wall *w, double *dw1, double *dw2);
 	// Differential equations describing chemical reactions inside the cells
 	virtual void CellDynamics(CellBase *c, double *dchem);
 	// to be executed after a cell division
 	virtual void OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2);
 	// to be executed for coloring a cell
 	virtual void SetCellColor(CellBase *c, QColor *color);
 	// return number of chemicals
 	virtual int NChem(void);
 };

src/TutorialCode/Tutorial3/mymodel.pro

➞

Show inline comments

@@ new file 100644 @@
+#
 # $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.
+#
 TARGET = mymodel
 VLEAFHOME = ../../..
 CONFIG += release
 CONFIG -= debug
 CONFIG += plugin
 BINDIR = $${VLEAFHOME}/bin
 LIBDIR = $${VLEAFHOME}/lib
 INCDIR = $${VLEAFHOME}/src
 DEFINES = QTGRAPHICS # VLEAFPLUGIN
 DESTDIR = $${BINDIR}/models
 HEADERS = $${TARGET}.h
 INCLUDEPATH += $${INCDIR}
 QMAKE_CXXFLAGS += -fexceptions #-I$${INCDIR}
 QMAKE_CXXFLAGS_DEBUG += -g3
 QMAKE_CXXFLAGS_DEBUG += -DQDEBUG
 QT += qt3support
 SOURCES = $${TARGET}.cpp
 TEMPLATE = lib
 unix {
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2
  QMAKE_LFLAGS += -fPIC
+}
 win32 {
  LIBXML2DIR = C:\libxml2
  LIBICONVDIR = C:\libiconv
  LIBZDIR = C:\libz
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -DLIBXML_STATIC
  QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include
+}
 # finish

src/TutorialCode/Tutorial4/mymodel.cpp

➞

Show inline comments

@@ new file 100644 @@
 /*
+ *
  *  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 <QObject>
 #include <QtGui>
 #include "simplugin.h"
 #include "parameter.h"
 #include "wallbase.h"
 #include "cellbase.h"
 #include "mymodel.h"
 static const std::string _module_id("$Id$");
 QString MyModel::ModelID(void) {
   // specify the name of your model here
   return QString( "Growth, division, coloring" );
+}
 // return the number of chemicals your model uses
 int MyModel::NChem(void) { return 0; }
 // To be executed after cell division
 void MyModel::OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2) {
   // rules to be executed after cell division go here
   // (e.g., cell differentiation rules)
+}
 void MyModel::SetCellColor(CellBase *c, QColor *color) {
   // add cell coloring rules here
 	if (c->Area()/c->BaseArea()>1.8) { color->setNamedColor("blue"); }
 	else { color->setNamedColor("green"); }
+}
 void MyModel::CellHouseKeeping(CellBase *c) {
   // add cell behavioral rules here
 	c->EnlargeTargetArea(par->cell_expansion_rate);
 	if (c->Area() > par->rel_cell_div_threshold * c->BaseArea()) {
 		c->Divide();
+	}
+}
 void MyModel::CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2) {
   // add biochemical transport rules here
+}
 void MyModel::WallDynamics(Wall *w, double *dw1, double *dw2) {
   // add biochemical networks for reactions occuring at walls here
+}
 void MyModel::CellDynamics(CellBase *c, double *dchem) {
   // add biochemical networks for intracellular reactions here
+}
 Q_EXPORT_PLUGIN2(mymodel, MyModel)

src/TutorialCode/Tutorial4/mymodel.h

➞

Show inline comments

@@ new file 100644 @@
 /*
  *  $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.
+ *
  */
 #include <QObject>
 #include <QtGui>
 #include <QString>
 #include "simplugin.h"
 class MyModel : public QObject, SimPluginInterface {
 	Q_OBJECT
 	Q_INTERFACES(SimPluginInterface);
 public:
 	virtual QString ModelID(void);
 	// Executed after the cellular mechanics steps have equillibrized
 	virtual void CellHouseKeeping (CellBase *c);
 	// Differential equations describing transport of chemicals from cell to cell
 	virtual void CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2);
 	// Differential equations describing chemical reactions taking place at or near the cell walls
 	// (e.g. PIN accumulation)
 	virtual void WallDynamics(Wall *w, double *dw1, double *dw2);
 	// Differential equations describing chemical reactions inside the cells
 	virtual void CellDynamics(CellBase *c, double *dchem);
 	// to be executed after a cell division
 	virtual void OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2);
 	// to be executed for coloring a cell
 	virtual void SetCellColor(CellBase *c, QColor *color);
 	// return number of chemicals
 	virtual int NChem(void);
 };

src/TutorialCode/Tutorial4/mymodel.pro

➞

Show inline comments

@@ new file 100644 @@
+#
 # $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.
+#
 TARGET = mymodel
 VLEAFHOME = ../../..
 CONFIG += release
 CONFIG -= debug
 CONFIG += plugin
 BINDIR = $${VLEAFHOME}/bin
 LIBDIR = $${VLEAFHOME}/lib
 INCDIR = $${VLEAFHOME}/src
 DEFINES = QTGRAPHICS # VLEAFPLUGIN
 DESTDIR = $${BINDIR}/models
 HEADERS = $${TARGET}.h
 INCLUDEPATH += $${INCDIR}
 QMAKE_CXXFLAGS += -fexceptions #-I$${INCDIR}
 QMAKE_CXXFLAGS_DEBUG += -g3
 QMAKE_CXXFLAGS_DEBUG += -DQDEBUG
 QT += qt3support
 SOURCES = $${TARGET}.cpp
 TEMPLATE = lib
 unix {
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2
  QMAKE_LFLAGS += -fPIC
+}
 win32 {
  LIBXML2DIR = C:\libxml2
  LIBICONVDIR = C:\libiconv
  LIBZDIR = C:\libz
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -DLIBXML_STATIC
  QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include
+}
 # finish

src/TutorialCode/Tutorial5/mymodel.cpp

➞

Show inline comments

@@ new file 100644 @@
 /*
+ *
  *  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 <QObject>
 #include <QtGui>
 #include "simplugin.h"
 #include "parameter.h"
 #include "wallbase.h"
 #include "cellbase.h"
 #include "mymodel.h"
 static const std::string _module_id("$Id$");
 QString MyModel::ModelID(void) {
   // specify the name of your model here
   return QString( "Growth hormones" );
+}
 // return the number of chemicals your model uses
 int MyModel::NChem(void) { return 1; }
 // To be executed after cell division
 void MyModel::OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2) {
   // rules to be executed after cell division go here
   // (e.g., cell differentiation rules)
 	// set one cell to source after first division
 	if (CellBase::NCells()==2) {
 		daughter1->SetCellType(1);
 		daughter2->SetCellType(0);
+	}
 	// if a source cells has divided, one of the daughters becomes the new source
 	if (daughter1->CellType()==1) {
 		// if both cells are at the tissue perimeter, choose at random
 		if (daughter1->AtBoundaryP() && daughter2->AtBoundaryP()) {
 			if (qrand()%2){
 				daughter1->SetCellType(1);
 				daughter2->SetCellType(0);
 			} else {
 				daughter1->SetCellType(0);
 				daughter2->SetCellType(1);
+			}
 		} else {
 			// otherwise choose the one that is still at the perimeter
 			if (daughter1->AtBoundaryP()) {
 				daughter1->SetCellType(1);
 				daughter2->SetCellType(0);
 			} else {
 				daughter1->SetCellType(0);
 				daughter2->SetCellType(1);
+			}
+		}
+	}
+}
 void MyModel::SetCellColor(CellBase *c, QColor *color) {
   // add cell coloring rules here
 	// white: high concentration of growth hormone, black low concentration
 	double val = c->Chemical(0)/(1.+c->Chemical(0));
 	color->setRgbF(val, val, val);
+}
 void MyModel::CellHouseKeeping(CellBase *c) {
   // add cell behavioral rules here
 	if (CellBase::NCells()==1)
 		// first cell expands unconditionally
 		c->EnlargeTargetArea(par->cell_expansion_rate);
 	else
 		c->EnlargeTargetArea(c->Chemical(0)*par->cell_expansion_rate);
 	if (c->Area() > par->rel_cell_div_threshold * c->BaseArea()) {
 		c->Divide();
+	}
+}
 void MyModel::CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2) {
   // add biochemical transport rules here
 	double phi = w->Length() * par->D[0] * ( w->C2()->Chemical(0) - w->C1()->Chemical(0) );
 	dchem_c1[0]+=phi;
 	dchem_c2[0]-=phi;
+}
 void MyModel::WallDynamics(Wall *w, double *dw1, double *dw2) {
   // add biochemical networks for reactions occuring at walls here
+}
 void MyModel::CellDynamics(CellBase *c, double *dchem) {
   // add biochemical networks for intracellular reactions here
 	if (c->CellType()==1) {
 		dchem[0] = par->leaf_tip_source;
+	}
+}
 Q_EXPORT_PLUGIN2(mymodel, MyModel)

src/TutorialCode/Tutorial5/mymodel.h

➞

Show inline comments

@@ new file 100644 @@
 /*
  *  $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.
+ *
  */
 #include <QObject>
 #include <QtGui>
 #include <QString>
 #include "simplugin.h"
 class MyModel : public QObject, SimPluginInterface {
 	Q_OBJECT
 	Q_INTERFACES(SimPluginInterface);
 public:
 	virtual QString ModelID(void);
 	// Executed after the cellular mechanics steps have equillibrized
 	virtual void CellHouseKeeping (CellBase *c);
 	// Differential equations describing transport of chemicals from cell to cell
 	virtual void CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2);
 	// Differential equations describing chemical reactions taking place at or near the cell walls
 	// (e.g. PIN accumulation)
 	virtual void WallDynamics(Wall *w, double *dw1, double *dw2);
 	// Differential equations describing chemical reactions inside the cells
 	virtual void CellDynamics(CellBase *c, double *dchem);
 	// to be executed after a cell division
 	virtual void OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2);
 	// to be executed for coloring a cell
 	virtual void SetCellColor(CellBase *c, QColor *color);
 	// return number of chemicals
 	virtual int NChem(void);
 };

src/TutorialCode/Tutorial5/mymodel.pro

➞

Show inline comments

@@ new file 100644 @@
+#
 # $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.
+#
 TARGET = mymodel
 VLEAFHOME = ../../..
 CONFIG += release
 CONFIG -= debug
 CONFIG += plugin
 BINDIR = $${VLEAFHOME}/bin
 LIBDIR = $${VLEAFHOME}/lib
 INCDIR = $${VLEAFHOME}/src
 DEFINES = QTGRAPHICS # VLEAFPLUGIN
 DESTDIR = $${BINDIR}/models
 HEADERS = $${TARGET}.h
 INCLUDEPATH += $${INCDIR}
 QMAKE_CXXFLAGS += -fexceptions #-I$${INCDIR}
 QMAKE_CXXFLAGS_DEBUG += -g3
 QMAKE_CXXFLAGS_DEBUG += -DQDEBUG
 QT += qt3support
 SOURCES = $${TARGET}.cpp
 TEMPLATE = lib
 unix {
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2
  QMAKE_LFLAGS += -fPIC
+}
 win32 {
  LIBXML2DIR = C:\libxml2
  LIBICONVDIR = C:\libiconv
  LIBZDIR = C:\libz
  LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -DLIBXML_STATIC
  QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include
+}
 # finish

src/VirtualLeafpar.tmpl

➞

Show inline comments

 title = Parameter values for The Virtual Leaf / title
 label = <b>Visualization</b> / label
 arrowcolor = white / string
 arrowsize = 100 / double
 textcolor = red / string
 cellnumsize = 1 / int
 nodenumsize = 1 / int
 node_mag = 1.0 / double
 outlinewidth = 1.0 / double
 cell_outline_color = forestgreen / string
-resize_stride = 10 / int
+resize_stride = 0 / int
 label = / label
 label = <b>Cell mechanics</b> / label
 T = 1.0 / double
 lambda_length = 100. / double
 lambda_celllength = 0. / double
 target_length = 60. / double
 cell_expansion_rate = 1. / double
 cell_div_expansion_rate = 0. / double
 auxin_dependent_growth = true / bool
 ode_accuracy = 1e-4 / double
 mc_stepsize = 0.4 / double
 mc_cell_stepsize = 0.2 / double
 energy_threshold = 1000. / double
 bend_lambda = 0. / double
 alignment_lambda = 0. / double
 rel_cell_div_threshold = 2. / double
 rel_perimeter_stiffness = 2 / double
 collapse_node_threshold = 0.05 / double
 morphogen_div_threshold = 0.2 / double
 morphogen_expansion_threshold = 0.01 / double
 copy_wall = true / bool
 label = / label
 label = <b>Auxin transport and PIN1 dynamics</b> / label
 source = 0. / double
 D = 0., 0.0, 0.0, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. / doublelist
 initval = 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. / doublelist
 k1 = 1. / double
 k2 = 0.3 / double
 r = 1. / double
 kr = 1. / double
 #k3 = 1 / double
 #k4 = 0.3 / double
 km = 1. / double
 Pi_tot = 1. / double
 transport = 0.036 / double
 ka = 1 / double
 # d = 1e-06 / double
 # e = 0.01 / double
 # f = 10 / double
 pin_prod = 0.001 / double
 pin_prod_in_epidermis = 0.1 / double
 pin_breakdown = 0.001 / double
 pin_breakdown_internal = 0.001 / double
 aux1prod = 0.001 / double
 aux1prodmeso = 0. / double
 aux1decay = 0.001 / double
 aux1decaymeso = 0.1 / double
 aux1transport = 0.036 / double
 aux_cons = 0. / double
 aux_breakdown = 0. / double
 kaux1 = 1 / double
 kap = 1 / double
 leaf_tip_source = 0.001 / double
 sam_efflux = 0.0001 / double
 sam_auxin = 10. / double
 # gf_prod = 1e-3 / double
 # gf_decay = 1e-3 / double
 sam_auxin_breakdown = 0 / double
 #label = / label
 #label = <b>miscellaneous</b> / label
 van3prod = 0.002 / double
 van3autokat = 0.1 / double
 van3sat = 10 / double
 k2van3 = 0.3 / double
 #glvprod = 0.0 / double
 #glvbreakdown = 0. / double
 label = / label
 label = <b>Integration parameters</b> / label
 dt = 0.1 / double
 rd_dt = 1.0 / double
 datadir = . / directory
 movie = false / bool
 nit = 100000 / int
 maxt = 1000. / double
 storage_stride = 10 / int
 xml_storage_stride = 500 / int
 rseed = -1 / int
 #label = / label
 #label = <b>Parameters for new chemical</b> / label
 #glvproduction = 0. / double
 #glvdecay = 0. / double
 #glvdiffthreshold = 0. / double
 label = / label
 label = <b>Meinhardt leaf venation model</b> / label
 constituous_expansion_limit = 16 / int
 vessel_inh_level = 1 / double
 vessel_expansion_rate = 0.25 / double
 d = 0. / double
 e = 0. / double
 f = 0. / double
 c = 0. / double
 mu = 0. / double
 nu = 0. / double
 rho0 = 0. / double
 rho1 = 0. / double
 c0 = 0. / double
 gamma = 0. / double
 eps = 0. / double
 label = / label
 label = <b>User-defined parameters</b> / label
 k = 0., 0.0, 0.0, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. / doublelist
 i1 = 0 / int
 i2 = 0 / int
 i3 = 0 / int
 i4 = 0 / int
 i5 = 0 / int
 s1 =  / string
 s2 =  / string
 s3 =  / string
 b1 = false / bool
 b2 = false / bool
 b3 = false / bool
 b4 = false / bool
 dir1 = . / directory
 dir2 = . / directory

src/build_models/auxingrowthplugin.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 <QObject>
 #include <QtGui>
 #include "../simplugin.h"
 #include "parameter.h"
 #include "wallbase.h"
 #include "cellbase.h"
 #include "auxingrowthplugin.h"
 #include "far_mem_5.h"
 static const std::string _module_id("$Id$");
 bool batch = false;
 // To be executed after cell division
-void AuxinGrowthPlugin::OnDivide(ParentInfo &parent_info, CellBase &daughter1, CellBase &daughter2) {
+void AuxinGrowthPlugin::OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2) {
 	// Auxin distributes between parent and daughter according to area
-	double area1 = daughter1.Area(), area2 = daughter2.Area();
+	double area1 = daughter1->Area(), area2 = daughter2->Area();
 	double tot_area = area1 + area2;
 	daughter1.SetChemical(0,daughter1.Chemical(0)*(area1/tot_area));
 	daughter2.SetChemical(0,daughter2.Chemical(0)*(area2/tot_area));
 	daughter1->SetChemical(0,daughter1->Chemical(0)*(area1/tot_area));
 	daughter2->SetChemical(0,daughter2->Chemical(0)*(area2/tot_area));
 	// After divisions, parent and daughter cells get a standard stock of PINs.
 	daughter1.SetChemical(1, par->initval[1]);
 	daughter2.SetChemical(1, par->initval[1]);
 	daughter1->SetChemical(1, par->initval[1]);
 	daughter2->SetChemical(1, par->initval[1]);
 	// Reset transporter values of parent and daughter
 	QList<WallBase *> walls;
 	foreach(WallBase *w, walls) {
-		w->setTransporter(&daughter1, 1, 0.);
 		w->setTransporter(daughter1, 1, 0.);
+	}
 	//daughter1.LoopWalls(Wall::setTransporter(&daughter1, 1, 0.));
 	/* for (list<Wall *>::const_iterator w=daughter2.walls.begin();
 		 w!=daughter2.walls.end();
 		 w++) {
 		// reset transporter value
 		(*w)->setTransporter(&daughter2, 1, 0.);
+	}
 	*/
+}
-void AuxinGrowthPlugin::SetCellColor(CellBase &c, QColor &color) {
+void AuxinGrowthPlugin::SetCellColor(CellBase *c, QColor *color) {
 	// Red: PIN1
 	// Green: Auxin
 	if (c.CellType()==1) color = QColor("Blue");
 	else color.setRgb(c.Chemical(1)/(1+c.Chemical(1)) * 255.,(c.Chemical(0)/(1+c.Chemical(0)) * 255.),/* (chem[2]/(1+chem[2]) *255.) */ 0);
 	if (c->CellType()==1) color->setNamedColor("Blue");
 	else color->setRgb(c->Chemical(1)/(1+c->Chemical(1)) * 255.,(c->Chemical(0)/(1+c->Chemical(0)) * 255.),/* (chem[2]/(1+chem[2]) *255.) */ 0);
+}
-void AuxinGrowthPlugin::CellHouseKeeping(CellBase &c) {
+void AuxinGrowthPlugin::CellHouseKeeping(CellBase *c) {
 	if (c.Boundary()==CellBase::None) {
 		if (c.Area() > par->rel_cell_div_threshold * c.BaseArea() ) {
 			c.SetChemical(0,0);
 			c.Divide();
 	if (c->Boundary()==CellBase::None) {
 		if (c->Area() > par->rel_cell_div_threshold * c->BaseArea() ) {
 			c->SetChemical(0,0);
 			c->Divide();
+        }
 		if (c.Chemical(0)>0.6) {
 			c.SetCellType(1);
 		if (c->Chemical(0)>0.6) {
 			c->SetCellType(1);
+		}
 		// expand according to auxin concentration
-   		c.EnlargeTargetArea(par->auxin_dependent_growth?(c.Chemical(0)/(1.+c.Chemical(0)))*par->cell_expansion_rate:par->cell_expansion_rate);
+   		c->EnlargeTargetArea(par->auxin_dependent_growth?(c->Chemical(0)/(1.+c->Chemical(0)))*par->cell_expansion_rate:par->cell_expansion_rate);
+	}
+}
 void AuxinGrowthPlugin::CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2) {
 	// leaf edge is const source of auxin
     // (Neumann boundary condition: we specify the influx)
     if (w->C2()->BoundaryPolP()) {
 		if (w->AuxinSource()) {
 			double aux_flux = par->leaf_tip_source * w->Length();
 			dchem_c1[0]+= aux_flux;
 			// dchem_c2 is undefined..!
 			return;
 		} else {
 			return;
+		}
+	}
 	if (w->C1()->BoundaryPolP()) {
 		if (w->AuxinSource()) {
 			double aux_flux = par->leaf_tip_source * w->Length();
 			dchem_c2[0] += aux_flux;
 			// dchem_c1 is undefined...!
 			return;
 		} else {
 			if (w->AuxinSink()) {
 				// efflux into Shoot Apical meristem
 				// we assume all PINs are directed towards shoot apical meristem
 				dchem_c2[0] -= par->sam_efflux * w->C2()->Chemical(0) / (par->ka + w->C2()->Chemical(0));
 				return;
 			} else
 				return;
+		}
+    }
 	// Passive fluxes (Fick's law)
     // only auxin flux now
     // flux depends on edge length and concentration difference
     int c=0;
     double phi = w->Length() * ( par->D[c] ) * ( w->C2()->Chemical(c) - w->C1()->Chemical(c) );
     dchem_c1[c] += phi;
     dchem_c2[c] -= phi;
     // Active fluxes (PIN1 mediated transport)
     // (Transporters measured in moles, here)
     // efflux from cell 1 to cell 2
     double trans12 = ( par->transport * w->Transporters1(1) * w->C1()->Chemical(0) / (par->ka + w->C1()->Chemical(0)) );
     // efflux from cell 2 to cell 1
     double trans21 = ( par->transport * w->Transporters2(1) * w->C2()->Chemical(0) / (par->ka + w->C2()->Chemical(0)) );
     dchem_c1[0] += trans21 - trans12;
     dchem_c2[0] += trans12 - trans21;
+}
 void AuxinGrowthPlugin::WallDynamics(Wall *w, double *dw1, double *dw2) {
     // Cells polarize available PIN1 to Shoot Apical Meristem
     if (w->C2()->BoundaryPolP()) {
 		if (w->AuxinSink()) {
 			dw1[0] = 0.; dw2[0] = 0.;
 			//dw1[2] = 0.; dw2[2] = 0.;
 			// assume high auxin concentration in SAM, to convince PIN1 to polarize to it
 			// exocytosis regulated0
 			double nb_auxin = par->sam_auxin;
 			double receptor_level = nb_auxin * par->r / (par->kr + nb_auxin);
 			dw1[1] = par->k1 * w->C1()->Chemical(1) * receptor_level /( par->km + w->C1()->Chemical(1) ) - par->k2 * w->Transporters1(1);
 			dw2[1] = 0.;
 			return;
 		} else {
 			dw1[0]=dw2[0]=dw1[1]=dw2[1];//=dw1[2]=dw2[2];
 			return;
+		}
+    }
     if (w->C1()->BoundaryPolP()) {
 		if (w->AuxinSink())  {
 			dw1[0] = 0.; dw2[0] = 0.;
 			//dw1[2] = 0.; dw2[2] = 0.;
 			// assume high auxin concentration in SAM, to convince PIN1 to polarize to it
 			// exocytosis regulated
 			double nb_auxin = par->sam_auxin;
 			double receptor_level = nb_auxin * par->r / (par->kr + nb_auxin);
 			dw2[1] = par->k1 * w->C2()->Chemical(1) * receptor_level /( par->km + w->C2()->Chemical(1) ) - par->k2 * w->Transporters2(1);
 			dw1[1] = 0.;
 			return;
 		}  else {
 			dw1[0]=dw2[0]=dw1[1]=dw2[1];//=dw1[2]=dw2[2];
 			return;
+		}
+    }
     // PIN1 localization at wall 1
     // Note: chemical 0 is Auxin (intracellular storage only)
     // Chemical 1 is PIN1 (walls and intracellular storage)
     //! \f$ \frac{d Pij/dt}{dt} = k_1 A_j \frac{P_i}{L_ij} - k_2 P_{ij} \f$
     // Note that Pij is measured in term of concentration (mol/L)
     // Pi in terms of quantity (mol)
     double dPijdt1=0., dPijdt2=0.;
     // normal cell
     double  auxin2 = w->C2()->Chemical(0);
     double receptor_level1 = auxin2 * par->r / (par->kr + auxin2);
     dPijdt1 =
 	// exocytosis regulated
     par->k1 * w->C1()->Chemical(1) * receptor_level1 / ( par->km + w->C1()->Chemical(1) ) - par->k2 * w->Transporters1(1);
     double  auxin1 = w->C1()->Chemical(0);
     double receptor_level2 = auxin1 * par->r / (par->kr + auxin1);
     // normal cell
     dPijdt2 =
 	// exocytosis regulated
 	par->k1 * w->C2()->Chemical(1) * receptor_level2 / ( par->km + w->C2()->Chemical(1) ) - par->k2 * w->Transporters2(1);
     /* PIN1 of neighboring vascular cell inhibits PIN1 endocytosis */
     dw1[0] = 0.; dw2[0] = 0.;
     //dw1[2] = 0.; dw2[2] = 0.;
     dw1[1] = dPijdt1;
     dw2[1] = dPijdt2;
+}
-double AuxinGrowthPlugin::complex_PijAj(CellBase &here, CellBase &nb, Wall &w) {
+double AuxinGrowthPlugin::complex_PijAj(CellBase *here, CellBase *nb, Wall *w) {
 	// gives the amount of complex "auxinreceptor-Pin1"  at the wall (at QSS)
 	//return here.Chemical(1) * nb.Chemical(0) / ( par->km + here.Chemical(1));
-	double nb_aux = (nb.BoundaryPolP() && w.AuxinSink()) ? par->sam_auxin : nb.Chemical(0);
+	double nb_aux = (nb->BoundaryPolP() && w->AuxinSink()) ? par->sam_auxin : nb->Chemical(0);
 	double receptor_level = nb_aux * par->r / (par->kr + nb_aux);
-	return here.Chemical(1) * receptor_level / ( par->km + here.Chemical(1));
+	return here->Chemical(1) * receptor_level / ( par->km + here->Chemical(1));
+}
 void AuxinGrowthPlugin::CellDynamics(CellBase *c, double *dchem) {
 	// Note: Pi and Pij measured in numbers of molecules, not concentrations
 		double dPidt = 0.;
 		double sum_Pij = c->SumTransporters( 1 );
 		// exocytosis regulated:
 	dPidt = -par->k1 * c->ReduceCellAndWalls<double>( far_3_arg_mem_fun( *this, &AuxinGrowthPlugin::complex_PijAj ) ) + par->k2 * sum_Pij;
 	/*for ( list<Wall *>::const_iterator w = c->walls.begin();
 		 w!=walls.end();
 		 w++) {
 		if ((*w)->C1() == c)
 			dPidt +=  complex_PijAj( (*w)->C1(), (*w)->C2(), *w );
 		else
 			dPidt +=  complex_PijAj( (*w)->C2(), (*w)->C1(), *w );
 	}*/
 	// production of PIN depends on auxin concentration
 	dPidt +=  (c->AtBoundaryP()?par->pin_prod_in_epidermis:par->pin_prod) * c->Chemical(0) - c->Chemical(1) * par->pin_breakdown;
 	// no PIN production in SAM
 	if (c->Boundary() == CellBase::SAM) {
 			dchem[1]=0.;
 		dchem[0]= - par->sam_auxin_breakdown * c->Chemical(0);
 	} else {
 		dchem[1] = dPidt;
 		// source of auxin
 		dchem[0] = par->aux_cons - par->aux_breakdown * c->Chemical(0);
+	}
+}
 Q_EXPORT_PLUGIN2(auxingrowthplugin, AuxinGrowthPlugin)

src/build_models/auxingrowthplugin.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.
+ *
  */
 #ifndef _AUXINGROWTHPLUGIN_H_
 #define _AUXINGROWTHPLUGIN_H_
 #include <QObject>
 #include <QtGui>
 #include <QString>
 #include "../simplugin.h"
 class AuxinGrowthPlugin : public QObject, SimPluginInterface {
 	Q_OBJECT
 	Q_INTERFACES(SimPluginInterface);
 public:
 	virtual QString ModelID(void) { return QString( "Auxin accumulation and growth" ); }
 	// Executed after the cellular mechanics steps have equillibrized
-	virtual void CellHouseKeeping (CellBase &c);
+	virtual void CellHouseKeeping (CellBase *c);
 	// Differential equations describing transport of chemicals from cell to cell
 	virtual void CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2);
 	// Differential equations describing chemical reactions taking place at or near the cell walls
 	// (e.g. PIN accumulation)
 	virtual void WallDynamics(Wall *w, double *dw1, double *dw2);
 	// Differential equations describing chemical reactions inside the cells
 	virtual void CellDynamics(CellBase *c, double *dchem);
 	// to be executed after a cell division
-	virtual void OnDivide(ParentInfo &parent_info, CellBase &daughter1, CellBase &daughter2);
+	virtual void OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2);
 	// to be executed for coloring a cell
-	virtual void SetCellColor(CellBase &c, QColor &color);
+	virtual void SetCellColor(CellBase *c, QColor *color);
 	// return number of chemicals
 	virtual int NChem(void) { return 2; }
 private:
-	double complex_PijAj(CellBase &here, CellBase &nb, Wall &w);
+	double complex_PijAj(CellBase *here, CellBase *nb, Wall *w);
 };
 #endif

src/build_models/meinhardtplugin.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 <QObject>
 #include <QtGui>
 #include "simplugin.h"
+#include "../simplugin.h"
 #include "parameter.h"
 #include "warning.h"
 #include "wallbase.h"
 #include "cellbase.h"
 #include "meinhardtplugin.h"
 static const std::string _module_id("$Id$");
 bool batch = false;
 // To be executed after cell division
-void MeinhardtPlugin::OnDivide(ParentInfo &parent_info, CellBase &daughter1, CellBase &daughter2) {
+void MeinhardtPlugin::OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2) {
+}
-void MeinhardtPlugin::SetCellColor(CellBase &c, QColor &color) {
+void MeinhardtPlugin::SetCellColor(CellBase *c, QColor *color) {
-	if (fpclassify(c.Chemical(0))==FP_NAN) {
+	if (fpclassify(c->Chemical(0))==FP_NAN) {
 		// somehow the function isnan doesn't work properly on my system... SuSE Linux
 		// 10.0 64-bits (isnan seems not be implemented using fpclassify).
 		MyWarning::warning("Whoops! Numerical instability!!");
-		color.setNamedColor("red");
+		color->setNamedColor("red");
 	} else {
 		double range_min = 0.;//, range_max = 1.;
-		if (c.Chemical(0)<range_min) {
+		if (c->Chemical(0)<range_min) {
 			MyWarning::warning("Whoops! Numerical instability!!");
-			color.setNamedColor("blue");
+			color->setNamedColor("blue");
 		} else {
-			color.setRgb(c.Chemical(1)/(1+c.Chemical(1)) * 255.,(c.Chemical(0)/(1+c.Chemical(0)) * 255.),(c.Chemical(3)/(1+c.Chemical(3)) *255.) );
+			color->setRgb(c->Chemical(1)/(1+c->Chemical(1)) * 255.,(c->Chemical(0)/(1+c->Chemical(0)) * 255.),(c->Chemical(3)/(1+c->Chemical(3)) *255.) );
+		}
+	}
+}
-void MeinhardtPlugin::CellHouseKeeping(CellBase &c) {
+void MeinhardtPlugin::CellHouseKeeping(CellBase *c) {
 	if (c.Area() > par->rel_cell_div_threshold * c.BaseArea() ) {
 		c.Divide();
 	if (c->Area() > par->rel_cell_div_threshold * c->BaseArea() ) {
 		c->Divide();
+	}
     // cell expansion is inhibited by substrate (chem 3)
 	if (!par->constituous_expansion_limit || c.NCells()<par->constituous_expansion_limit) {
 		c.EnlargeTargetArea(par->cell_expansion_rate );
 	if (!par->constituous_expansion_limit || c->NCells()<par->constituous_expansion_limit) {
 		c->EnlargeTargetArea(par->cell_expansion_rate );
 	} else {
-		if (c.Chemical(0)<0.5) {
+		if (c->Chemical(0)<0.5) {
 			double tmp;
-			c.EnlargeTargetArea((tmp=(1.-par->vessel_inh_level*c.Chemical(3))*par->cell_expansion_rate /* + c.Chemical(4)*/)<0?0:tmp);
+			c->EnlargeTargetArea((tmp=(1.-par->vessel_inh_level*c->Chemical(3))*par->cell_expansion_rate /* + c->Chemical(4)*/)<0?0:tmp);
 		} else {
-			c.EnlargeTargetArea(par->vessel_expansion_rate);
+			c->EnlargeTargetArea(par->vessel_expansion_rate);
+		}
+	}
+}
 void MeinhardtPlugin::CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2) {
 	// No flux boundaries for all chemicals, except activator: boundary is sink
   if (w->C1()->BoundaryPolP() || w->C2()->BoundaryPolP()) {
 		if (w->C1()->BoundaryPolP()) {
 			dchem_c2[1] -=  w->Length() * ( par->D[1] ) * ( w->C2()->Chemical(1) );
 			} else {
 			dchem_c1[1] -=  w->Length() * ( par->D[1] ) * ( w->C1()->Chemical(1) );
+			}
 		return;
+		}
     // Passive fluxes (Fick's law)
     for (int c=0;c<NChem();c++) {
 		double phi = w->Length() * ( par->D[c] ) * ( w->C2()->Chemical(c) - w->C1()->Chemical(c) );
 		dchem_c1[c] += phi;
 		dchem_c2[c] -= phi;
+	}
+}
 void MeinhardtPlugin::WallDynamics(Wall *w, double *dw1, double *dw2) {
 	for (int c = 0;c<NChem();c++) {
 		dw1[c] = 0.; dw2[c] = 0.;
+    }
+}
 void MeinhardtPlugin::CellDynamics(CellBase *c, double *dchem) {
 	double Y = c->Chemical(0);
 	double A = c->Chemical(1);
 	double H = c->Chemical(2);
 	double S = c->Chemical(3);
 	//double expansin = c->Chemical(4);
     dchem[0] = ( par->d * A - par->e * Y + Y*Y/(1 + par->f * Y*Y ) );
 	dchem[1] = ( par->c * A*A*S/H - par->mu * A + par->rho0*Y );
 	dchem[2] = ( par->c * A*A*S - par->nu*H + par->rho1*Y );
 	dchem[3] = ( par->c0 - par->gamma*S - par->eps * Y * S );
 	//dchem[4] = ( -par->expansindecay * expansin );
 	//for (int i=0;i<4;i++) { cerr << "[ " << dchem[i] << " ]"; } cerr << endl;
 //	cerr << "Chemicals: "; for (int i=0;i<NChem();i++) { cerr << c->Chemical(i) << " "; } cerr << endl;
 // test:	dchem[0] = 0.01 * c->Chemical(0) * ( 1. - c->Chemical(0));
+}
 Q_EXPORT_PLUGIN2(meinhardtplugin, MeinhardtPlugin)

src/build_models/meinhardtplugin.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.
+ *
  */
 #ifndef _MEINHARDTPLUGIN_H_
 #define _MEINHARDTPLUGIN_H_
 #include <QObject>
 #include <QtGui>
 #include <QString>
 #include "simplugin.h"
+#include "../simplugin.h"
 class MeinhardtPlugin : public QObject, SimPluginInterface {
 	Q_OBJECT
 	Q_INTERFACES(SimPluginInterface);
 public:
 	virtual QString ModelID(void) { return QString( "Meinhardt 1976, with growth" ); }
 	// Executed after the cellular mechanics steps have equillibrized
-	virtual void CellHouseKeeping (CellBase &c);
+	virtual void CellHouseKeeping (CellBase *c);
 	// Differential equations describing transport of chemicals from cell to cell
 	virtual void CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2);
 	// Differential equations describing chemical reactions taking place at or near the cell walls
 	// (e.g. PIN accumulation)
 	virtual void WallDynamics(Wall *w, double *dw1, double *dw2);
 	// Differential equations describing chemical reactions inside the cells
 	virtual void CellDynamics(CellBase *c, double *dchem);
 	// to be executed after a cell division
-	virtual void OnDivide(ParentInfo &parent_info, CellBase &daughter1, CellBase &daughter2);
+	virtual void OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2);
 	// to be executed for coloring a cell
-	virtual void SetCellColor(CellBase &c, QColor &color);
+	virtual void SetCellColor(CellBase *c, QColor *color);
 	// return number of chemicals
 	virtual int NChem(void) { return 4; }
 };
 #endif

src/build_models/plugin_auxingrowth.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
 LIBDIR = ../../lib
 DEFINES = QTGRAPHICS # VLEAFPLUGIN
 DESTDIR = $${BINDIR}/models
 TARGET = auxingrowth
 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
+ LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2
  QMAKE_LFLAGS += -fPIC
+}
 win32 {
  LIBXML2DIR = C:\libxml2
  LIBICONVDIR = C:\libiconv
  LIBZDIR = C:\libz
  LIBS += -Llib -lvleaf
+ LIBS += -L$${LIBDIR} -Llib -lvleaf
  QMAKE_CXXFLAGS += -DLIBXML_STATIC
  QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include
+}
 # finis

src/build_models/plugin_leaf.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
 LIBDIR = ../../lib
 DEFINES = QTGRAPHICS # VLEAFPLUGIN
 DESTDIR = $${BINDIR}/models
 TARGET = leaf
 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
+ LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2
  QMAKE_LFLAGS += -fPIC
+}
 win32 {
  LIBXML2DIR = C:\libxml2
  LIBICONVDIR = C:\libiconv
  LIBZDIR = C:\libz
  LIBS += -Llib -lvleaf
+ LIBS += -L$${LIBDIR} -Llib -lvleaf
  QMAKE_CXXFLAGS += -DLIBXML_STATIC
  QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include
+}
 # finis

src/build_models/plugin_meinhardt.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
 LIBDIR = ../../lib
 DEFINES = QTGRAPHICS # VLEAFPLUGIN
 DESTDIR = $${BINDIR}/models
 TARGET = meinhardt
 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
+ LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2
  QMAKE_LFLAGS += -fPIC
+}
 win32 {
  LIBXML2DIR = C:\libxml2
  LIBICONVDIR = C:\libiconv
  LIBZDIR = C:\libz
  LIBS += -Llib -lvleaf
+ LIBS += -L$${LIBDIR} -Llib -lvleaf
  QMAKE_CXXFLAGS += -DLIBXML_STATIC
  QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include
+}
 # finis

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
 LIBDIR = ../../lib
 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
+ LIBS += -L$${LIBDIR} -lvleaf
  QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2
  QMAKE_LFLAGS += -fPIC
+}
 win32 {
  LIBXML2DIR = C:\libxml2
  LIBICONVDIR = C:\libiconv
  LIBZDIR = C:\libz
  LIBS += -Llib -lvleaf
+ LIBS += -L$${LIBDIR} -Llib -lvleaf
  QMAKE_CXXFLAGS += -DLIBXML_STATIC
  QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include
+}
 # finis

src/build_models/translate_plugin.pl

➞

Show inline comments

@@ new file 100755 @@
 #!/usr/bin/perl
 $cfilename = shift(@ARGV) || die "Usage: translate_plugin.pl [cfile] [hfile] [profile]\n";
 $hfilename = shift(@ARGV) || die "Usage: translate_plugin.pl [cfile] [hfile] [profile]\n";
 $pfilename = shift(@ARGV) || die "Usage: translate_plugin.pl [cfile] [hfile] [profile]\n";
 $ocfname = $cfilename; $ocfname =~ s/\.cpp/_tl.cpp/g;
 $ohfname = $hfilename; $ohfname =~ s/\.h/_tl.h/g;
 $opfname = $pfilename; $opfname =~ s/\.pro/_tl.pro/g;
 print STDERR "Translating '$cfilename' to '$ocfname', '$hfilename' to '$ohfname', and '$pfilename' to '$opfname'\n";
 open cfile,"<$cfilename";
 open ocfile,">$ocfname";
 while (<cfile>) {
 	#s/$hfilename/$ohfname/g;
 	# translate function definitions
 	if (/[a-zA-Z0-9 ]*::OnDivide/) {
 		s/ParentInfo &parent_info/ParentInfo *parent_info/g;
 		s/CellBase &daughter1/CellBase *daughter1/g;
 		s/CellBase &daughter2/CellBase *daughter2/g;
+	}
 	if (/[a-zA-Z0-9 ]*::SetCellColor/) {
 		s/CellBase &c/CellBase *c/g;
 		s/QColor &color/QColor *color/g;
+	}
 	if (/[a-zA-Z0-9 ]*::CellHouseKeeping/) {
 		s/CellBase &c/CellBase *c/g;
+	}
 	# translate member function calls
 	s/\bparent_info\b\./parent_info->/g;
 	s/\bdaughter1\b\./daughter1->/g;
 	s/\bdaughter2\b\./daughter2->/g;
 	s/\bc\b\./c->/g;
 	s/\bcolor\b\./color->/g;
 	print ocfile;
+}
 open hfile,"<$hfilename";
 open ohfile,">$ohfname";
 while (<hfile>) {
 	if (/[ \t]*virtual[ \t]+void[ \t]+CellHouseKeeping/) {
 		s/CellBase &c/CellBase *c/g;
+	}
 	if (/[ \t]*virtual[ \t]+void[ \t]+OnDivide/) {
 		s/ParentInfo &parent_info/ParentInfo *parent_info/g;
 		s/CellBase &daughter1/CellBase *daughter1/g;
 		s/CellBase &daughter2/CellBase *daughter2/g;
+	}
 	if (/[ \t]*virtual[ \t]+void[ \t]+SetCellColor/) {
 		s/CellBase &c/CellBase *c/g;
 		s/QColor &color/QColor *color/g;
+	}
 	print ohfile;
+}
 open pfile,"<$pfilename";
 open opfile,">$opfname";
 while (<pfile>) {
 	s/\bplugin\b\.h/plugin_tl\.h/g;
 	s/\bplugin\b\.cpp/plugin_tl\.cpp/g;
 	print opfile;
+}
@@ \ No newline at end of file @@

src/cell.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 <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];
+	}
 	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();
 			if ((div_edges[i].first->BoundaryP() && div_edges[i].second->BoundaryP()) && // Both edge nodes are boundary nodes AND
 			     ((m->findNextBoundaryNode(div_edges[i].first))->Index() == div_edges[i].second->Index())){ // The boundary proceeds from first to second.
                                 #ifdef QDEBUG
 			        qDebug() << "Index of the first node: " << div_edges[i].first->Index() << endl;
 			        qDebug() << "Index of the second node: " << div_edges[i].second->Index() << endl;
 			        qDebug() << "Boundary proceeds from: " <<  div_edges[i].first->Index()
 				         << "to: " << (m->findNextBoundaryNode(div_edges[i].first))->Index() << endl << endl;
                                 #endif
 			        new_node_ind[i]->SetBoundary();
 				// We will need to repair the boundary polygon later, since we will insert new nodes
 				//cerr << "Boundary touched for Node " << new_node_ind[i]->Index() << "\n";
 				boundary_touched_flag=true;
 				// and insert it into the boundary_polygon
 				// find the position of the first node in the boundary
 				list<Node *>::iterator ins_pos = find
 				(m->boundary_polygon->nodes.begin(),
 				 m->boundary_polygon->nodes.end(),
 				 div_edges[i].first);
 				// ... second node comes before or after it ...
 				if (*(++ins_pos!=m->boundary_polygon->nodes.end()?
 					  ins_pos:m->boundary_polygon->nodes.begin())!=div_edges[i].second) {
 					m->boundary_polygon->nodes.insert(((ins_pos--)!=m->boundary_polygon->nodes.begin()?ins_pos:(--m->boundary_polygon->nodes.end())), new_node_ind[i]);
 					// .. set the neighbors of the new node ...
 					// in this case e.second and e.first are inverted
 				} else {
 					// insert before second node, so leave ins_pos as it is,
 					// that is: incremented
 					m->boundary_polygon->nodes.insert(ins_pos, new_node_ind[i]);
 					// .. set the neighbors of the new node ...
+				}
+			}
 			list<Neighbor> owners;
 			// push all cells owning the two nodes of the divides edges
 			// onto a list
 			copy((div_edges[i].first)->owners.begin(),
 				 (div_edges[i].first)->owners.end(),
 				 back_inserter(owners));
 			copy((div_edges[i].second)->owners.begin(),
 				 (div_edges[i].second)->owners.end(),
 				 back_inserter(owners));
 			// find first non-self duplicate in the owners:
 			// cells owning the same two nodes
 			// share an edge with me
 			owners.sort( mem_fun_ref( &Neighbor::Cmp ) );
                         #ifdef QDEBUG
 			list<Neighbor> unique_owners;
 			copy(owners.begin(), owners.end(), back_inserter(unique_owners));
 			unique_owners.unique( mem_fun_ref( &Neighbor::Eq ) );
 			qDebug() << "The dividing edge nodes: " << div_edges[i].first->Index() << " and " << div_edges[i].second->Index() << " are owned by cells: ";
 			// spit out each owners' cell index
 			foreach(Neighbor neighbor, unique_owners){
 			  qDebug() << neighbor.cell->Index() << "  ";
+			}
 			qDebug() << endl;
                         #endif
 			// Search through the sorted list of edge node owners looking for duplicate pairs. Each pair represents an actual edge owner.
 			list<Neighbor> edge_owners;
 			list<Neighbor>::iterator it;
 			for (it=owners.begin(); it!=owners.end(); it++) {
 			  it = adjacent_find(it, owners.end(), neighbor_cell_eq);
 			  if (it == owners.end()) break; // bail if reach the end of the list
                           #ifdef QDEBUG
 			  qDebug() << "Considering: " << it->cell->Index() << " as a possible edge owner." << endl;
                           #endif
 			  if (it->cell->Index() != this->Index()) {
                             #ifdef QDEBUG
 			    qDebug() << "Adding: " << it->cell->Index() << " to the list of edge owners." << endl;
                             #endif
 			    edge_owners.push_back(*it);
+			  }
+			}
 			if (edge_owners.size() > 1){
 			  // Remove the boundary polygon - if its there
 			  list<Neighbor>::iterator it;
 			  if ((it = find_if (edge_owners.begin(), edge_owners.end(), bind2nd(mem_fun_ref(&Neighbor::CellEquals), -1))) != edge_owners.end()) {
                             #ifdef QDEBUG
 			    qDebug() << "deleating: " << it->cell->Index() << " from the list of edge owners." << endl;
                             #endif
 			    edge_owners.erase(it);
+			  }
+			}
                         #ifdef QDEBUG
 			qDebug() << "The edge owners list has: " << edge_owners.size() << " elements" << endl;
 			#endif
 			// Since the list should always contain exactly one element, pass it on as an iterator
 			list<Neighbor>::iterator c = (edge_owners.size() != 0) ? edge_owners.begin() : edge_owners.end();
 			// (can we have more than one neighboring cell here??)
 			if (c!=owners.end()) {
 				neighbor_cell = c->cell;
 				if (!c->cell->BoundaryPolP()) {
 					// find correct position in the cells node list
 					// to insert the new node
 					list<Node *>::iterator ins_pos = find
 					(neighbor_cell->nodes.begin(),
 					 neighbor_cell->nodes.end(),
 					 div_edges[i].first);
 					neighbor_cell->nodes.insert(ins_pos, new_node_ind[i]);
 					neighbor_cell->ConstructConnections();
 					// give walls to daughter later
+			  }
 			} else {
 				neighbor_cell = 0;
+			}
+		}
 		// Split the Wall with the neighboring cell
 		// if the neighbor cell has not yet been identified above, do it now
 		if (neighbor_cell == 0) {
 			list<Neighbor> owners;
 			// push all cells owning the two nodes of the divides edges
 			// onto a list
 			copy((div_edges[i].first)->owners.begin(),
 				 (div_edges[i].first)->owners.end(),
 				 back_inserter(owners));
 			copy((div_edges[i].second)->owners.begin(),
 				 (div_edges[i].second)->owners.end(),
 				 back_inserter(owners));
 			// find first non-self duplicate in the owners:
 			// cells owning the same two nodes
 			// share an edge with me
 			owners.sort( mem_fun_ref( &Neighbor::Cmp ) );
 			list<Neighbor>::iterator c;
 			for (c=owners.begin();
 				 c!=owners.end();
 				 c++) {
 				c=adjacent_find(c,owners.end(),neighbor_cell_eq);
 				if (c->cell->Index() != this->Index() || c==owners.end()) break;
+			}
 			if (c!=owners.end())
 				neighbor_cell = c->cell;
 			else
 				neighbor_cell = 0;
+		}
 		if (neighbor_cell /* && !neighbor_cell->BoundaryPolP() */) {
 			//cerr << "Cell "  << index << " says: neighboring cell is " << neighbor_cell->index << endl;
 			/*************** 1. Find the correct wall element  ********************/
 			list<Wall *>::iterator w, start_search;
 			w = start_search = walls.begin();
 			do {
 				// Find wall between this cell and neighbor cell
 				w = find_if( start_search, walls.end(), bind2nd (mem_fun( &Wall::is_wall_of_cell_p ), neighbor_cell ) );
 				start_search = w; start_search++; // continue searching at next element
 			} while ( w!=walls.end() && !(*w)->IntersectsWithDivisionPlaneP( from, to ) ); // go on until we find the right one.
 			if (w == walls.end()) {
 				cerr << "Whoops, wall element not found...!\n";
 				cerr << "Cell ID: " << neighbor_cell->Index() << endl;
 				cerr << "My cell ID: " << Index() << endl;
 			} else {
 				// 2. Split it up, if we should (sometimes, the new node coincides with an existing node so
 				// we should not split up the Wall)
 				if (new_node_ind[i]!=(*w)->n1 && new_node_ind[i]!=(*w)->n2) {
 					Wall *new_wall;
 					// keep the length of the original wall; we need it to equally divide the transporter concentrations
 					// over the two daughter walls
 					(*w)->SetLength(); // make sure we've got the current length
 					orig_length[i] = (*w)->Length();
 					//cerr << "Original length is " << orig_length[i] << endl;
 					if ((*w)->c1 == this ) {
 						//  cerr << "Cell " << (*w)->c1->Index() << " splits up wall " << *(*w) << ", into: " << endl;
 						new_wall = new Wall( (*w)->n1, new_node_ind[i], this, neighbor_cell);
 						(*w)->n1 = new_node_ind[i];
 						//  cerr << "wall " << *(*w) << ", and new wall " << *new_wall << endl;
 					} else {
 						new_wall = new Wall( (*w)->n1, new_node_ind[i], neighbor_cell, this);
 						(*w)->n1 = new_node_ind[i];
+					}
 					//new_wall->ResetTransporterConcentrations(orig_length);
 					//(*w)->ResetTransporterConcentrations(orig_length);
 					// reset the transporter concentrations
 					/*	  new_wall->SetLength();
 					 new_wall->CorrectLength(orig_length);
 					 (*w)->SetLength();
 					 (*w)->CorrectLength(orig_length);*/
 					// 3. Give wall elements to appropriate cells
 					if (new_wall->n1 != new_wall->n2) {
 						if (par.copy_wall)
 							new_wall->CopyWallContents(**w);
 						else {
 							// If wall contents are not copied, decide randomly which wall will be the "parent"
 							// otherwise we will get biases (to the left), for example in the meristem growth model
 							if (RANDOM()<0.5) {
 								new_wall->SwapWallContents(*w);
+							}
+						}
 						AddWall(new_wall);
 						// cerr << "Building new wall: this=" << Index() << ", neighbor_cell = " << neighbor_cell->Index() << endl;
 						neighbor_cell->AddWall( new_wall);
 						//cerr << "Existing wall: c1 = " << (*w)->c1->Index() << ", neighbor_cell = " << (*w)->c2->Index() << endl;
 						// Remember the addresses of the new walls
 						div_wall[2*i+0] = *w;
 						div_wall[2*i+1] = new_wall;
 						// we will correct the transporter concentrations later in this member function, after division
 						// First the new nodes should be inserted into the cells, before we can calculate wall lengths
 						// Remember that cell walls can be bent, so have a bigger length than the Euclidean distance n1->n2
 					} else {
 						delete new_wall;
+					}
+				}
+			}
+		}
 	}  // closing loop over the two divided edges (for (int i=0;i<2;i++) )
 	// move half of the nodes to the daughter
+	{
 		//cerr << "Daughter: ";
 		list<Node *>::iterator start, stop;
 		start=new_node_locations.front();
 		//cerr << "*new_node_locations.front() = " << *new_node_locations.front() << endl;
 		if (new_node_flag[0]==1) {
 			start++;
 			if (start==nodes.end())
 				start=nodes.begin();
+		}
 		stop=new_node_locations.back();
 		if (new_node_flag[1]==2) {
 			if (stop==nodes.begin())
 				stop=nodes.end();
 			stop--;
+		}
 		list<Node *>::iterator i=start;
 		while ( i!=stop) {
 			// give the node to the daughter
 			// (find references to parent cell from this node,
 			// and remove them)
 			list<Neighbor>::iterator neighb_with_this_cell=
 			find_if((*i)->owners.begin(),
 					(*i)->owners.end(),
 				bind2nd(mem_fun_ref( &Neighbor::CellEquals ),this->Index() )  );
 			if (neighb_with_this_cell==(*i)->owners.end()) {
 				cerr << "not found\n";
 				abort();
+			}
 			(*i)->owners.erase(neighb_with_this_cell);
 			daughter->nodes.push_back( *i );
 			i++;
 			if (i==nodes.end())
 				i=nodes.begin();
 		};
+	}
 	// new node list of parent
 	list<Node *> new_nodes_parent;
 	// half of the nodes stay with the parent
+	{
 		list<Node *>::iterator start, stop;
 		start=new_node_locations.back();
 		if (new_node_flag[1]==1) {
 			start++;
 			if (start==nodes.end())
 				start=nodes.begin();
+		}
 		stop=new_node_locations.front();
 		if (new_node_flag[0]==2) {
 			if (stop==nodes.begin())
 				stop=nodes.end();
 			stop--;
+		}
 		list<Node *>::iterator i=start;
 		while (i!=stop) {
 			new_nodes_parent.push_back( *i );
 			i++;
 			if (i==nodes.end())
 				i = nodes.begin();
 		};
+	}
 	// insert shared wall
 	// insert shared nodes on surface of parent cell
 	new_nodes_parent.push_back( new_node_ind[0] );
 	daughter->nodes.push_back ( new_node_ind[1] );
 	// optionally add the new node to the nodeset (passed by pointer)
 	// (in this way we can move the NodeSet as a whole; useful for a fixed cutting line)
 	if (node_set) {
 		node_set->AddNode( new_node_ind[0] );
+	}
 #define MULTIPLE_NODES
 #ifdef MULTIPLE_NODES
 	// intermediate, extra nodes
 	// Calculate distance between the two new nodes
 	double dist=( new_node[1] - new_node[0] ).Norm();
 	//bool fixed_wall = (new_node_ind[0])->fixed && (new_node_ind[1])->fixed;
 	bool fixed_wall = false;
 	// Estimate number of extra nodes in wall
 	// factor 4 is to keep tension on the walls;
 	// this is a hidden parameter and should be made explicit
 	// later on.
 	int n=(int)((dist/Node::target_length)/4+0.5);
 	Vector nodevec = ( new_node[1]- new_node[0]).Normalised();
 	double element_length = dist/(double)(n+1);
 	// note that wall nodes need to run in inverse order in parent
 	list<Node *>::iterator ins_pos = daughter->nodes.end();
 	for (int i=1;i<=n;i++) {
 		Node *node=
 		m->AddNode( new Node( new_node[0] + i*element_length*nodevec ) );
 		node->fixed=fixed_wall;
 		if (!fix_cellwall)
 			node->boundary = false;
 		else { // if fix_cellwall is true, that is if we are cutting off
 			// part of a leaf to make a nice initial condition, we also want to make it part of the boundary
 			//node->boundary = true;
 			node->fixed = true;
 			//node->sam = true;
+		}
 		ins_pos=daughter->nodes.insert(ins_pos, node );
 		new_nodes_parent.push_back( node );
 		// optionally add the new node to the nodeset (passed by pointer)
 		// (in this way we can move the NodeSet as a whole; useful for a fixed cutting line)
 		if (node_set) {
 			node_set->AddNode( node );
+		}
+	}
 #endif
 	daughter->nodes.push_back( new_node_ind[0] );
 	new_nodes_parent.push_back( new_node_ind[1] );
 	// optionally add the new node to the nodeset (passed by pointer)
 	// (in this way we can move the NodeSet as a whole; useful for a fixed cutting line)
 	if (node_set) {
 		node_set->AddNode( new_node_ind[1] );
+	}
 	// move the new nodes to the parent
 	nodes.clear();
 	copy( new_nodes_parent.begin(),
 		 new_nodes_parent.end(),
 		 back_inserter(nodes) );
 	// Repair cell lists of Nodes, and node connectivities
 	ConstructConnections();
 	daughter->ConstructConnections();
 	if (boundary_touched_flag) {
 		m->boundary_polygon->ConstructConnections();
+	}
 	// collecting neighbors of divided cell
 	list<CellBase *> broken_neighbors;
 	// this cell's old neighbors
 	copy(neighbors.begin(), neighbors.end(), back_inserter(broken_neighbors) );
 	// this cell
 	broken_neighbors.push_back(this);
 	// its daughter
 	broken_neighbors.push_back(daughter);
 	// Recalculate area of parent and daughter
 	area = CalcArea();
 	daughter->area = daughter->CalcArea();
 	SetIntegrals();
 	daughter->SetIntegrals();
 	// Add a "Cell Wall" for diffusion algorithms
 	Wall *wall = new Wall( new_node_ind[0], new_node_ind[1], this, daughter );
 	AddWall( wall );
 	daughter->AddWall( wall );
 	//cerr << "Correct walls of cell " << Index() << " and daughter " << daughter->Index() << endl;
 	// Move Walls to daughter cell
 	list <Wall *> copy_walls = walls;
 	for (list<Wall *>::iterator w = copy_walls.begin();
 		 w!=copy_walls.end();
 		 w++) {
 		//cerr << "Doing wall, before:  " << **w << endl;
 		//  checks the nodes of the wall and gives it away if appropriate
 		(*w)->CorrectWall ( );
 		//cerr << "and after: " << **w << endl;
+	}
 	// Correct tranporterconcentrations of divided walls
 	for (int i=0;i<4;i++) {
 		if (div_wall[i]) {
 			div_wall[i]->SetLength();
 			div_wall[i]->CorrectTransporters(orig_length[i/2]);
+		}
+	}
 	//neighbors.push_back( daughter );
 	//daughter->neighbors.push_back( this );
 	//cerr << "Cell " << index << " has been dividing, and gave birth to Cell " << daughter->index << endl;
 	// now reconstruct neighbor list for all "broken" neighbors
 	for (list<CellBase *>::iterator i=broken_neighbors.begin();
 		 i!=broken_neighbors.end();i++) {
 		((Cell *)(*i))->ConstructNeighborList();
+	}
 	ConstructNeighborList();
 	daughter->ConstructNeighborList();
-	m->plugin->OnDivide(parent_info,*daughter, *this);
+	m->plugin->OnDivide(&parent_info, daughter, this);
 	// wall->OnWallInsert();
 	//daughter->OnDivide();
 	daughter->div_counter=(++div_counter);
+}
 // Move the whole cell
 void Cell::Move(const Vector T) {
     for (list<Node *>::const_iterator i=nodes.begin();
 		 i!=nodes.end();
 		 i++) {
 		*(*i)+=T;
+    }
+}
 double Cell::Displace(double dx, double dy, double dh) {
 	// Displace whole cell, add resulting energy to dh,
 	// and accept displacement if energetically favorable
 	//
 	// Method is called if a "fixed" node is displaced
 	// Warning: length constraint not yet  CORRECTLY implemented for this function
 	// Attempt to move this cell in a random direction
 	//  Vector movement(par.mc_cell_stepsize*(RANDOM()-0.5),par.mc_cell_stepsize*(RANDOM()-0.5),0);
 	dh=0;
 	Vector movement(dx,dy,0);
 	vector< pair<Node *, Node *> > length_edges;
 	vector<double> cellareas;
 	cellareas.reserve(neighbors.size());
 	// for the length constraint, collect all edges to this cell's nodes,
 	// which are not part of the cell
 	// the length of these edges will change
 	double old_length=0.;
 	for (list<Node *>::const_iterator i=nodes.begin();
 		 i!=nodes.end();
 		 i++) {
 		//if ((*i)->Fixed()) return; // commented out 01/12/05
 		for (list<Neighbor>::const_iterator n=(*i)->owners.begin();
 			 n!=(*i)->owners.end();
 			 n++) {
 			if (n->getCell()!=this) {
 				//if (!(m->getNode(n->nb1).Fixed() && m->getNode(n->nb2).Fixed())) {
 				length_edges.push_back( pair <Node *,Node *> (*i, n->nb1) );
 				length_edges.push_back( pair <Node *,Node *> (*i, n->nb2) );
 				old_length +=
 				DSQR(Node::target_length-(*(*i)-*(n->nb1)).Norm())+
 				DSQR(Node::target_length-(*(*i)-*(n->nb2)).Norm());
 				//}
+			}
+		}
+	}
 	// calculate area energy difference of neighboring cells
 	// (this cells' shape remains unchanged)
 	double old_area_energy=0., old_length_energy=0.;
 	for (list<CellBase *>::const_iterator i=neighbors.begin();
 		 i!=neighbors.end();
 		 i++) {
 		old_area_energy += DSQR((*i)->Area()-(*i)->TargetArea());
 		old_length_energy += DSQR((*i)->Length()-(*i)->TargetLength());
+	}
 	Move(movement);
 	double new_area_energy=0., new_length_energy=0.;
 	for (list<CellBase *>::const_iterator i=neighbors.begin();
 		 i!=neighbors.end();
 		 i++) {
 		cellareas.push_back((*i)->CalcArea());
 		new_area_energy += DSQR(cellareas.back()-(*i)->TargetArea());
 		new_length_energy += DSQR((*i)->CalcLength()-(*i)->TargetLength());
+	}
 	double new_length=0;
 	for ( vector< pair< Node *, Node * > >::const_iterator e = length_edges.begin();
 		 e != length_edges.end();
 		 e++) {
 		new_length +=  DSQR(Node::target_length-
 							(*(e->first)-*(e->second)).Norm());
+	}
 	dh += (new_area_energy - old_area_energy) + (new_length_energy - old_length_energy) * lambda_celllength +
 	par.lambda_length * (new_length - old_length);
 	if (dh<0 || RANDOM()<exp(-dh/par.T)) {
 		// update areas of cells
 		//cerr << "neighbors: ";
 		list<CellBase *>::const_iterator nb_it = neighbors.begin();
 		for (vector<double>::const_iterator ar_it = cellareas.begin();
 			 ar_it!=cellareas.end();
 			 ( ar_it++, nb_it++) ) {
 			((Cell *)(*nb_it))->area = *ar_it;
 			(*nb_it)->SetIntegrals();
+		}
 		//cerr << endl;
 		/*vector<double> area1;
 		 vector<double> area2;
 		 m->ExtractFromCells( mem_fun_ref(&Cell::Area), back_inserter(area1) );
 		 m->ExtractFromCells( mem_fun_ref(&Cell::CalcArea), back_inserter(area2));
 		 vector<double>::iterator i=area1.begin();
 		 vector<double>::iterator j=area2.begin();
 		 int c=0;
 		 for (;
 		 i!=area1.end();
 		 (i++, j++)) {
 		 if ( (*i-*j) > 1e-10) {
 		 cerr << c++ << " " << *i << " " << *j << endl;
 		 abort();
+		 }
 		 }*/
 	} else {
 		Move ( -1*movement);
+	}
 	return dh;
+}
 void Cell::Displace (void) {
 	Displace(par.mc_cell_stepsize*(RANDOM()-0.5),par.mc_cell_stepsize*(RANDOM()-0.5),0);
+}
 // Get energy level of whole cell (excluding length constraint?)
 double Cell::Energy(void) const {
 	double energy = 0.;
 	double length_contribution = 0.;
 	for (list<Node *>::const_iterator i=nodes.begin();
 		 i!=nodes.end();
 		 i++) {
 		for (list<Neighbor>::const_iterator n=(*i)->owners.begin();
 			 n!=(*i)->owners.end();
 			 n++) {
 			if (n->getCell()==this) {
 				length_contribution +=
 				DSQR(Node::target_length-(*(*i)-*(n->nb1)).Norm())+
 				DSQR(Node::target_length-(*(*i)-*(n->nb2)).Norm());
+			}
+		}
+	}
 	// wall elasticity constraint
 	energy += par.lambda_length * length_contribution;
 	// area constraint
 	energy += DSQR(CalcArea() - target_area);
 	// cell length constraint
 	energy += lambda_celllength * DSQR(Length() - target_length);
 	return energy;
+}
 bool Cell::SelfIntersect(void) {
     // The (obvious) O(N*N) algorithm
     // Compare each edge against each other edge
     // An O(N log(N)) algorithm by Shamos & Hoey (1976) supposedly exists;
     // it was mentioned on comp.graphics.algorithms
     // But I haven't been able to lay my hand on the paper.
     // Let's try whether we need it....
     // method used: http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
     for (list<Node *>::const_iterator i=nodes.begin();
 		 i!=nodes.end();
 		 i++) {
 		list<Node *>::const_iterator j=i;
 		++j;
 		for (;
 			 j!=nodes.end();
 			 j++)
+		{
 			Vector v1 = *(*i);
 			list<Node *>::const_iterator nb=i;
 			nb++;
 			if (nb == nodes.end()) {
 				nb = nodes.begin();
+			}
 			Vector v2 = *(*nb);
 			Vector v3 = *(*j);
 			nb=j;
 			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);*/
 			if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) {
 				//cerr << "ua = " << ua << ", ub = " << ub << endl;
 				return true;
+			}
+		}
+    }
     return false;
+}
 bool Cell::MoveSelfIntersectsP(Node *moving_node_ind, Vector new_pos) {
 	// Check whether the polygon will self-intersect if moving_node_ind
 	// were displaced to new_pos
 	// Compare the two new edges against each other edge
 	// O(2*N)
 	// method used for segment intersection:
 	// http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
 	Vector neighbor_of_moving_node[2];
 	//cerr << "list<Node *>::const_iterator moving_node_ind_pos = find (nodes.begin(),nodes.end(),moving_node_ind);\n";
 	list<Node *>::const_iterator moving_node_ind_pos = find (nodes.begin(),nodes.end(),moving_node_ind);
 	list<Node *>::const_iterator nb = moving_node_ind_pos;
 	//cerr << "Done\n";
 	nb++;
 	if (nb == nodes.end()) {
 		nb = nodes.begin();
+	}
 	neighbor_of_moving_node[0]=*(*nb);
 	nb=moving_node_ind_pos;
 	if (nb == nodes.begin()) {
 		nb = nodes.end();
+	}
 	nb--;
 	neighbor_of_moving_node[1]=*( *nb );
 	for (list<Node *>::const_iterator i=nodes.begin();
 		 i!=nodes.end();
 		 i++) {
 		for (int j=0;j<2;j++) { // loop over the two neighbors of moving node
 			list<Node *>::const_iterator nb=i;
 			nb++;
 			if (nb == nodes.end()) {
 				nb = nodes.begin();
+			}
 			if (*i == moving_node_ind || *nb == moving_node_ind) {
 				// do not compare to self
 				continue;
+			}
 			Vector v3 = *(*i);
 			Vector v4 = *(*nb);
 			double denominator =
 			(v4.y - v3.y)*(neighbor_of_moving_node[j].x - new_pos.x) - (v4.x - v3.x)*(neighbor_of_moving_node[j].y - new_pos.y);
 			double ua =
 			((v4.x - v3.x)*(new_pos.y - v3.y) - (v4.y - v3.y)*(new_pos.x -v3.x))/denominator;
 			double ub =
 			((neighbor_of_moving_node[j].x - new_pos.x)*(new_pos.y-v3.y) - (neighbor_of_moving_node[j].y- new_pos.y)*(new_pos.x - v3.x))/denominator;
 			/* double intersec_x = new_pos.x + ua*(neighbor_of_moving_node[j].x-new_pos.x);
 			 double intersec_y = new_pos.y + ua*(neighbor_of_moving_node[j].y-new_pos.y);*/
 			if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) {
 				//cerr << "ua = " << ua << ", ub = " << ub << endl;
 				return true;
+			}
+		}
+	}
 	return false;
+}
 /*! \brief Test if this cell intersects with the given line.
  */
 bool Cell::IntersectsWithLineP(const Vector v1, const Vector v2) {
 	// Compare the line against each edge
 	// method used: http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
 	for (list<Node *>::const_iterator i=nodes.begin();
 		 i!=nodes.end();
 		 i++)
+    {
 		Vector v3 = *(*i);
 		list<Node *>::const_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);*/
 		if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) {
 			return true;
+		}
+    }
 	return false;
+}
 /*! \brief Constructs Walls, but only one per cell boundary.
  Standard method constructs a Wall for each cell wall element,
  making transport algorithms computationally more intensive than needed.
  We can remove this? Well, let's leave it in the code in case we need it for something else. E.g. for importing leaf architectures in different formats than our own... :-)
  */
 void Cell::ConstructWalls(void) {
 	return;
 	if (dead) return;
 	walls.clear();
 	neighbors.clear();
 	// Get "corner points; i.e. nodes where more than 2 cells are connected
 	list<Node *> corner_points;
 	for (list<Node *>::const_iterator i=nodes.begin();
 		 i!=nodes.end();i++) {
 		// look for nodes belonging to >2 cells
 		if ((*i)->owners.size()>2) {
 			// push onto list
 			corner_points.push_back(*i);
+		}
+	}
 	// Construct Walls between corner points
 	// previous one in list
 	list<Node *>::const_iterator nb = (--corner_points.end());
 	// loop over list,
 	for (list<Node *>::const_iterator i=corner_points.begin();
 		 i!=corner_points.end(); ( i++, nb++) ) {
 		if (nb==corner_points.end()) nb=corner_points.begin();
 		// add owning cells to a list
 		list<Cell *> owning_cells;
 		Node &n(*(*i));
 		for (list<Neighbor>::const_iterator j=n.owners.begin();
 			 j!=n.owners.end();
 			 j++) {
 			owning_cells.push_back(j->cell);
+		}
 		Node &n2(*(*nb));
 		for (list<Neighbor>::const_iterator j=n2.owners.begin();
 			 j!=n2.owners.end();
 			 j++) {
 			owning_cells.push_back(j->cell);
+		}
 		// sort cell owners
 		owning_cells.sort( mem_fun( &Cell::Cmp ));
 		// find duplicates
 		vector<Cell *> duplicates;
 		list<Cell *>::const_iterator prevj = (--owning_cells.end());
 		for (list<Cell *>::const_iterator j=owning_cells.begin();
 			 j!=owning_cells.end();
 			 ( j++, prevj++) ) {
 			if (prevj==owning_cells.end()) prevj=owning_cells.begin();
 			if (*j==*prevj) duplicates.push_back(*j);
+		}
 		if (duplicates.size()==3) { // ignore cell boundary (this occurs only after the first division, I think)
 			vector<Cell *>::iterator dup_it = find_if(duplicates.begin(),duplicates.end(),mem_fun(&Cell::BoundaryPolP) );
 			if (dup_it!=duplicates.end())
 				duplicates.erase(dup_it);
 			else {
 				return;
+			}
+		}
 		// One Wall for each neighbor, so we should be able to correctly construct neighbor lists here.
 		if (duplicates[0]==this) {
 			//walls. new Wall(*nb,*i,duplicates[0],duplicates[1]) );
 			AddWall(  new Wall(*nb,*i,duplicates[0],duplicates[1]) );
 			if (!duplicates[1]->BoundaryPolP()) {
 				neighbors.push_back(duplicates[1]);
+			}
 		} else {
 			//walls.push_back( new Wall(*nb,*i,duplicates[1],duplicates[0]) );
 			AddWall ( new Wall(*nb,*i,duplicates[1],duplicates[0]) );
 			if (!duplicates[0]->BoundaryPolP()) {
 				neighbors.push_back(duplicates[0]);
+			}
+		}
+	}
+}
 void BoundaryPolygon::Draw(QGraphicsScene *c, QString tooltip) {
 	// Draw the BoundaryPolygon on a QCanvas object
 	CellItem* p = new CellItem(this, c);
 	QPolygonF pa(nodes.size());
 	int cc=0;
 	for (list<Node *>::const_iterator n=nodes.begin();
 		 n!=nodes.end();
 		 n++) {
 		Node *i=*n;
 		pa[cc++] = QPoint((int)((Offset().x+i->x)*Factor()),
 						  (int)((Offset().y+i->y)*Factor()) );
+	}
 	p->setPolygon(pa);
 	p->setPen(par.outlinewidth>=0?QPen( QColor(par.cell_outline_color),par.outlinewidth):QPen(Qt::NoPen));
 	p->setBrush( Qt::NoBrush );
 	p->setZValue(1);
 	if (!tooltip.isEmpty())
 		p->setToolTip(tooltip);
 	p->show();
+}
 void Cell::Flux(double *flux, double *D)  {
 	// loop over cell edges
 	for (int c=0;c<NChem();c++) flux[c]=0.;
 	for (list<Wall *>::iterator i=walls.begin();
 		 i!=walls.end();
 		 i++) {
 		// leaf cannot take up chemicals from environment ("no flux boundary")
 		if ((*i)->c2->BoundaryPolP()) continue;
 		// flux depends on edge length and concentration difference
 		for (int c=0;c<NChem();c++) {
 			double phi = (*i)->length * ( D[c] ) * ( ((Cell *)(*i)->c2)->chem[c] - chem[c] );
 			if ((*i)->c1!=this) {
 				cerr << "Warning, bad cells boundary: " << (*i)->c1->Index() << ", " << index << endl;
+			}
 			flux[c] += phi;
+		}
+	}
+}
 // graphics stuff, not compiled for batch versions
 #ifdef QTGRAPHICS
 #include "canvas.h"
 void Cell::Draw(QGraphicsScene *c, QString tooltip) {
 	// Draw the cell on a QCanvas object
 	if (DeadP()) {
 		cerr << "Cell " << index << " not drawn, because dead.\n";
 		return;
+	}
 	CellItem* p = new CellItem(this, c);
 	QPolygonF pa(nodes.size());
 	int cc=0;
 	for (list<Node *>::const_iterator n=nodes.begin();
 		 n!=nodes.end();
 		 n++) {
 		Node *i=*n;
 		pa[cc++] = QPoint((int)((offset[0]+i->x)*factor),
 						  (int)((offset[1]+i->y)*factor) );
+	}
 	QColor cell_color;
-	m->plugin->SetCellColor(*this,cell_color);
+	m->plugin->SetCellColor(this,&cell_color);
 	p->setPolygon(pa);
 	p->setPen(par.outlinewidth>=0?QPen( QColor(par.cell_outline_color),par.outlinewidth):QPen(Qt::NoPen));
 	p->setBrush( cell_color );
 	p->setZValue(1);
 	if (!tooltip.isEmpty())
 		p->setToolTip(tooltip);
 	p->show();
+}
 void Cell::DrawCenter(QGraphicsScene *c) const {
   // Maginfication derived similarly to that in nodeitem.cpp
   // Why not use Cell::Magnification()?
   const double mag = par.node_mag;
 	// construct an ellipse
   QGraphicsEllipseItem *disk = new QGraphicsEllipseItem ( -1*mag, -1*mag, 2*mag, 2*mag, 0, c );
 	disk->setBrush( QColor("forest green") );
 	disk->setZValue(5);
 	disk->show();
 	Vector centroid=Centroid();
 	disk -> setPos((offset[0]+centroid.x)*factor,(offset[1]+centroid.y)*factor);
+}
 void Cell::DrawNodes(QGraphicsScene *c) const {
 	for (list<Node *>::const_iterator n=nodes.begin();
 		 n!=nodes.end();
 		 n++) {
 		Node *i=*n;
 		//QCanvasEllipse *item = new QCanvasEllipse( 10, 10, c);
 		NodeItem *item = new NodeItem ( &(*i), c );
 		//QGraphicsRectItem *item = new QGraphicsRectItem(-50, -50, 50, 50, 0, c);
 		//disk->setBrush( QColor("IndianRed") );
 		/*if (i->sam) {
 		 item->setBrush( purple );
 		 } else {
 		 if (i->boundary) {
 		 item->setBrush( deep_sky_blue );
+		 }
 		 else {
 		 item->setBrush( indian_red );
+		 }
 		 }*/
 		item->setColor();
 		/*(if (item->getNode().DeadP()) {
 		 item->setBrush( QBrush (Qt::Dense6Pattern) );
 		 }*/
 		item->setZValue(5);
 		item->show();
 		item ->setPos(((offset[0]+i->x)*factor),
 					  ((offset[1]+i->y)*factor) );
+	}
+}
 void Cell::DrawIndex(QGraphicsScene *c) const {
 	//  stringstream text;
 	//     text << index;
 	//     Vector centroid = Centroid();
 	//     QCanvasText *number = new QCanvasText ( QString (text.str()), c );
 	//     number->setColor( QColor(par.textcolor) );
 	//     number->setZ(20);
 	//     number->setFont( QFont( "Helvetica", par.cellnumsize, QFont::Bold) );
 	//     number->show();
 	//     number -> move((int)((offset[0]+centroid.x)*factor),
 	// 		   (int)((offset[1]+centroid.y)*factor) );
 	DrawText( c, QString("%1").arg(index));
+}
 // Draw any text in the cell's center
 void Cell::DrawText(QGraphicsScene *c, const QString &text) const {
 	Vector centroid = Centroid();
 	QGraphicsSimpleTextItem *ctext = new QGraphicsSimpleTextItem ( text, 0, c );
 	ctext->setPen( QPen(QColor(par.textcolor)) );
 	ctext->setZValue(20);
 	ctext->setFont( QFont( "Helvetica", par.cellnumsize, QFont::Bold) );
 	//ctext->setTextFlags(Qt::AlignCenter);
 	ctext->show();
 	ctext ->setPos(((offset[0]+centroid.x)*factor),
 				   ((offset[1]+centroid.y)*factor) );
+}
 void Cell::DrawAxis(QGraphicsScene *c) const {
 	Vector long_axis;
 	double width;
 	Length(&long_axis, &width);
 	//cerr << "Length is "  << length << endl;
 	long_axis.Normalise();
 	Vector short_axis=long_axis.Perp2D();
 	Vector centroid = Centroid();
 	Vector from = centroid - 0.5 * width * short_axis;
 	Vector to = centroid + 0.5 * width *short_axis;
 	QGraphicsLineItem *line = new QGraphicsLineItem(0, c);
 	line->setPen( QPen(QColor(par.arrowcolor),2) );
 	line->setZValue(2);
 	line->setLine( ( (offset[0]+from.x)*factor ),
 				  ( (offset[1]+from.y)*factor ),
 				  ( (offset[0]+to.x)*factor ),
 				  ( (offset[1]+to.y)*factor ) );
 	line->setZValue(10);
 	line->show();
+}
 void Cell::DrawStrain(QGraphicsScene *c) const {
 	MyWarning::warning("Sorry, Cell::DrawStrain temporarily not implemented.");
 	/* Vector long_axis;
 	double width;
 	Length(&long_axis, &width);
 	//cerr << "Length is "  << length << endl;
 	long_axis.Normalise();
 	Vector short_axis=long_axis.Perp2D();
 	//  To test method "Strain" temporarily substitute "short_axis" for "strain"
 	Vector strain = Strain();
 	//strain.Normalise();
 	//static ofstream strainf("strain.dat");
 	//strainf << strain.Norm() << endl;
 	Vector centroid = Centroid();
 	// Vector from = centroid - 0.5 * width * short_axis;
     // Vector to = centroid + 0.5 * width *short_axis;
 	Vector from = centroid - 0.5 * strain;
 	Vector to = centroid + 0.5 * strain;
 	QGraphicsArrowItem *arrow = new QGraphicsArrowItem(0, c);
 	arrow->setPen( QPen(QColor(par.arrowcolor),100) );
 	arrow->setLine( ( (offset[0]+from.x)*factor ),
 				   ( (offset[1]+from.y)*factor ),
 				   ( (offset[0]+to.x)*factor ),
 				   ( (offset[1]+to.y)*factor ) );
 	arrow->setZValue(10.);
 	arrow->show();
 	*/
+}
 // Draw connecting lines to neighbors
 /*void Cell::DrawTriangles(QCanvas &c) {
  for (list<Neighbor>::const_iterator nb=nb_list.begin();
  nb!=nb_list.end();
  nb++) {
  QCanvasLine *line = new QCanvasLine(&c);
  line->setPen( QPen(QColor("black"),2) );
  line->setZ(2);
  line->setPoints((offset[0]+x)*factor,(offset[1]+y)*factor,
  (offset[0]+nb->c->x)*factor,(offset[1]+nb->c->y)*factor);
  line->setZ(10);
  line->show();
+ }
  }*/
 void Cell::DrawFluxes(QGraphicsScene *c, double arrowsize)  {
 	// get the mean flux through this cell
 	//Vector vec_flux = ReduceWalls( mem_fun_ref( &Wall::VizFlux ), Vector() );
 	Vector vec_flux = ReduceCellAndWalls<Vector>( PINdir );
 	vec_flux.Normalise();
 	vec_flux *= arrowsize;
 	QGraphicsArrowItem *arrow = new QGraphicsArrowItem(0,c);
 	Vector centroid = Centroid();
 	Vector from = centroid - vec_flux/2.;
 	Vector to = centroid + vec_flux/2.;
 	arrow->setPen( QPen(QColor(par.arrowcolor),par.outlinewidth));
 	arrow->setZValue(2);
 	arrow->setLine( ( (offset[0]+from.x)*factor ),
 				   ( (offset[1]+from.y)*factor ),
 				   ( (offset[0]+to.x)*factor ),
 				   ( (offset[1]+to.y)*factor ) );
 	arrow->setZValue(10);
 	arrow->show();
+}
 void Cell::DrawWalls(QGraphicsScene *c) const {
 	for_each(walls.begin(), walls.end(), bind2nd ( mem_fun ( &Wall::Draw ) , c ) );
 	// to see the cells connected the each wall (for debugging), uncomment the following
 	//for_each(walls.begin(), walls.end(), bind2nd ( mem_fun ( &Wall::ShowStructure ), c ) );
+}
 void Cell::DrawValence(QGraphicsScene *c) const {
 	DrawText(c, QString("%1").arg(walls.size()) );
+}
 #endif
 /*! \brief Recalculate the lengths of the cell's Walls.
  Call this function after the Monte Carlo updates, and before doing the reaction-diffusion iterations.
  */
 void Cell::SetWallLengths(void) {
 	for (list<Wall *>::iterator de=walls.begin();
 		 de!=walls.end();
 		 de++) {
 		// Step 1: find the path of nodes leading along the Wall.
 		// A Wall often represents a curved cell wall: we want the total
 		// length _along_ the wall here...
 		// Locate first and second nodes of the edge in list of nodes
 		list<Node *>::const_iterator first_node_edge = find(nodes.begin(), nodes.end(), (*de)->n1);
 		list<Node *>::const_iterator second_node_edge_plus_1 = ++find(nodes.begin(), nodes.end(), (*de)->n2);
 		double sum_length = 0.;
 		// Now, walk to the second node of the edge in the list of nodes
 		for (list<Node *>::const_iterator n=++first_node_edge;
 			 n!=second_node_edge_plus_1;
 			 ++n ) {
 			if (n==nodes.end()) n=nodes.begin(); /* wrap around */
 			list<Node *>::const_iterator prev_n = n;
 			if (prev_n==nodes.begin()) prev_n=nodes.end();
 			--prev_n;
 			// Note that Node derives from a Vector, so we can do vector calculus as defined in vector.h
 			sum_length += (*(*prev_n) - *(*n)).Norm();
 			//cerr << "Node " << *prev_n << " to " << *n << ", cumulative length = " << sum_length << endl;
+		}
 		// We got the total length of the Wall now, store it:
 		(*de)->length = sum_length;
 		//cerr << endl;
 		// goto next de
+	}
+}
 //! Add Wall w to the list of Walls
 void Cell::AddWall( Wall *w ) {
 	// if necessary, we could try later inserting it at the correct position
 	if (w->c1 == w->c2 ){
 		cerr << "Wall between identical cells: " << w->c1->Index()<< endl;
+	}
 	// Add Wall to Cell's list
 	walls.push_back( w );
 	// Add wall to Mesh's list if it isn't there yet
 	if (find (
 			  m->walls.begin(), m->walls.end(),
 			  w )
 		== m->walls.end() ) {
 		m->walls.push_back(w);
+	}
+}
 //! Remove Wall w from the list of Walls
 list<Wall *>::iterator Cell::RemoveWall( Wall *w ) {
 	// remove wall from Mesh's list
 	m->walls.erase(
 				   find(
 						m->walls.begin(), m->walls.end(),
 						w )
 				   );
 	// remove wall from Cell's list
 	return
 	walls.erase (
 				 find(
 					  walls.begin(), walls.end(),
 					  w )
 				 );
+}
 void Cell::EmitValues(double t) {
 	//  cerr << "Attempting to emit " << t << ", " << chem[0] << ", " << chem[1] << endl;
 	//chem[3] = SumTransporters( 1 );
 	emit ChemMonValue(t, chem);
+}

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
 	if (!lambda_celllength) {
 		// Without length constraint we do not keep track of the cells'
 		// moments of inertia. So we must calculate them here.
 		SetIntegrals();
+	}
 	double intrx=intgrl_x/6.;
 	double intry=intgrl_y/6.;
 	double ixx=(intgrl_xx/12.)-(intrx*intrx)/area;
 	double ixy=(intgrl_xy/24.)+(intrx*intry)/area;
 	double iyy=(intgrl_yy/12.)-(intry*intry)/area;
 	double rhs1=(ixx+iyy)/2., rhs2=sqrt( (ixx-iyy)*(ixx-iyy)+4*ixy*ixy )/2.;
 	double lambda_b=rhs1+rhs2;
 	// see: http://scienceworld.wolfram.com/physics/MomentofInertiaEllipse.html
 	//    cerr << "n = " << n << "\n";
 	// Vector eigenvectors[2];
 	// eigenvectors[0] = Vector(-(-ixx + iyy ) + rhs2, ixy, 0);
 	// eigenvectors[1] = Vector(-(-ixx + iyy ) - rhs2, ixy, 0);
 	if (long_axis) {
 		*long_axis = Vector(-ixy, lambda_b - ixx, 0);
 		//   cerr << "ixx = " << ixx << ", ixy = " << ixy << ", iyy = " << iyy << ", area = " << area << endl;
+	}
 	if (width) {
 		*width = 4*sqrt((rhs1-rhs2)/area);
+	}
 	return 4*sqrt(lambda_b/area);
+}
 double CellBase::CalcLength(Vector *long_axis, double *width)  const {
 	// Calculate length and axes of CellBase, without touching cells raw moments
 	// Calculate inertia tensor
 	// see file inertiatensor.nb for explanation of this method
 	double my_intgrl_xx=0., my_intgrl_xy=0., my_intgrl_yy=0.;
 	double my_intgrl_x=0., my_intgrl_y=0., my_area=0.;
 	my_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();
 		my_area+=(*i)->x*(*nb)->y;
 		my_area-=(*nb)->x*(*i)->y;
 		my_intgrl_xx+=
 			((*i)->x*(*i)->x+
 			 (*nb)->x*(*i)->x+
 			 (*nb)->x*(*nb)->x ) *
 			((*i)->x*(*nb)->y-
 			 (*nb)->x*(*i)->y);
 		my_intgrl_xy+=
 			((*nb)->x*(*i)->y-
 			 (*i)->x*(*nb)->y)*
 			((*i)->x*(2*(*i)->y+(*nb)->y)+
 			 (*nb)->x*((*i)->y+2*(*nb)->y));
 		my_intgrl_yy+=
 			((*i)->x*(*nb)->y-
 			 (*nb)->x*(*i)->y)*
 			((*i)->y*(*i)->y+
 			 (*nb)->y*(*i)->y+
 			 (*nb)->y*(*nb)->y );
 		my_intgrl_x+=
 			((*nb)->x+(*i)->x)*
 			((*i)->x*(*nb)->y-
 			 (*nb)->x*(*i)->y);
 		my_intgrl_y+=
 			((*nb)->y+(*i)->y)*
 			((*i)->x*(*nb)->y-
 			 (*nb)->x*(*i)->y);
+	}
 	//my_area/=2.0;
 	my_area = fabs(my_area)/2.0;
 	double intrx=my_intgrl_x/6.;
 	double intry=my_intgrl_y/6.;
 	double ixx=(my_intgrl_xx/12.)-(intrx*intrx)/my_area;
 	double ixy=(my_intgrl_xy/24.)+(intrx*intry)/my_area;
 	double iyy=(my_intgrl_yy/12.)-(intry*intry)/my_area;
 	double rhs1=(ixx+iyy)/2., rhs2=sqrt( (ixx-iyy)*(ixx-iyy)+4*ixy*ixy )/2.;
 	double lambda_b=rhs1+rhs2;
 	// see: http://scienceworld.wolfram.com/physics/MomentofInertiaEllipse.html
 	//    cerr << "n = " << n << "\n";
 	// Vector eigenvectors[2];
 	// eigenvectors[0] = Vector(-(-ixx + iyy ) + rhs2, ixy, 0);
 	// eigenvectors[1] = Vector(-(-ixx + iyy ) - rhs2, ixy, 0);
 	if (long_axis) {
 		*long_axis = Vector(-ixy, lambda_b - ixx, 0);
 		//   cerr << "ixx = " << ixx << ", ixy = " << ixy << ", iyy = " << iyy << ", my_area = " << my_area << endl;
+	}
 	if (width) {
 		*width = 4*sqrt((rhs1-rhs2)/my_area);
+	}
 	return 4*sqrt(lambda_b/my_area);
+}
 // void CellBase::NodeRemoved(int n) {
 //   for (list<Node *>::iterator i=nodes.begin();
 //        i!=nodes.end();
 //        i++) {
 //     if ((*i)->Index()>n) {
 //       (*i)->index--;
 //     }
 //   }
 // }
 void CellBase::ConstructNeighborList(void) {
 	neighbors.clear();
 	for (//list<Wall *>::const_reverse_iterator wit=walls.rbegin();
 		 list<Wall *>::const_iterator wit=walls.begin();
 		 // somehow the reverse_iterator returns by walls needs to be casted to const to let this work.
 		 // it seems to me it is a bug in the STL implementation...
 		 //wit!=(list<Wall *>::const_reverse_iterator)walls.rend();
 		 wit!=walls.end();
 		 wit++) {
 		if ((*wit)->C1() != this) {
 			neighbors.push_back((*wit)->C1());
 		} else {
 			neighbors.push_back((*wit)->C2());
+		}
+	}
 	/*
 	 for (list<CellBase *>::iterator e=neighbors.begin();
 		  e!=neighbors.end();
 		  e++) {
 		 cerr << (*e)->Index() << " ";
 		 if ((*e)->CellBase::BoundaryPolP()) {
 			 cerr << " b ";
+		 }
+	 }
 	 */
 	// remove all boundary_polygons from the list
 	list <CellBase *>::iterator e=neighbors.begin();
 	at_boundary=false;
 	do {
 		// Code crashes here after cutting off part of the leaf. I can't find the problem.
 		// Leaving the "Illegal" walls in the simulation helps. (c1=-1 && c2=-1)
 		// Work-around: define leaf primordium. Save to XML. Restart. Read XML file.
 		// Sorry about this; I hope to solve this annoying issue later. RM :-).
 		// All cells in neighbors seem to be okay (I might be messing some part of the memory elsewhere
 		// during the cutting operation?).
 		e = find_if(neighbors.begin(),neighbors.end(),mem_fun(&CellBase::BoundaryPolP));
 		if (e!=neighbors.end()) {
 			e=neighbors.erase(e);
 			at_boundary=true;
 		} else {
 			break;
+		}
   } while(1);
+}
 // CellBase constructs its neighbor list from its node lists
 // Assumes, obviously, that the node lists are up to date
 // (i.e. call ConstructConnections before calling this method)
 // We'll keep this one private, anyway.
 /* void CellBase::ConstructNeighborList(void) {
 //  extern ofstream debug_stream;
 neighbors.clear();
 //  debug_stream << "Nodes: ";
 //  copy(nodes.begin(),nodes.end(),ostream_iterator<Node>(debug_stream, " "));
 //debug_stream << endl;
 for (list<Node *>::const_iterator i=nodes.begin();
 	 i!=nodes.end();
 	 i++) {
     // collect all cells to which my nodes are connected on one list
     //transform((*i)->cells.begin(),(*i)->cells.end(), back_inserter(neighbors), mem_fun_ref(&Neighbor::CellBase));
     // index of next node
     list<Node *>::const_iterator nn=i;
     ++nn;
     if (nn==nodes.end())
 		nn=nodes.begin();
     //    debug_stream << "Node " << *i << ", Neighbor " << *nn << endl;
     // debug_stream << "Owners: ";
     //    copy((*i)->cells.begin(),(*i)->cells.end(),ostream_iterator<Neighbor>(debug_stream, " "));
     // debug_stream << endl;
     for (list<Neighbor>::const_iterator nb=(*i)->owners.begin();
 		 nb!=(*i)->owners.end();
 		 nb++) {
 		// collect info about neighboring cells, not about myself
 		if (nb->CellBase!=this) {
 			// make sure the whole edge touches this putative neighbor
 			// if (*nn == nb->nb1 || *nn == nb->nb2) {
 			//walls.push_back( new Wall(*i,*nn,this,nb->CellBase) );
 			//debug_stream << "Adding edge " << walls.back() << " to CellBase " << index << endl;
 			//}
 			neighbors.push_back( nb->CellBase );
+		}
+    }
+}
 neighbors.sort();
 list<CellBase *>::iterator e=unique(neighbors.begin(),neighbors.end());
 // iterator e point to the end of the subsequence of unique elements
 // remove every thing that comes after it
 neighbors.erase(e, neighbors.end());
 // okay, now neighbors contains all neighbors of this CellBase, including itself
 // A future optimization for the diffusion algorithm: now we list
 // each of the edges of a (curved) CellBase boundary separately.  We
 // could keep track just of the total length per CellBase boundary
 // the following is not necessary anymore. Is
 // checked at earlier stage
 // // remove myself from the list
 //   e = find(neighbors.begin(),neighbors.end(),index);
 //   if (e!=neighbors.end())
 //   neighbors.erase(e);
 //
 // remove boundary_polygon from the list (CellBase identity <0 )
 e=neighbors.begin();
 at_boundary=false;
 do {
     e = find_if(neighbors.begin(),neighbors.end(),mem_fun(&CellBase::BoundaryPolP));
     if (e!=neighbors.end()) {
 		e=neighbors.erase(e);
 		at_boundary=true;
     } else {
 		break;
+    }
 } while(1);
 }*/
 /*void Cell::print_nblist(void) const {
 //  cerr << "{ ";
 for (list<Neighbor>::const_iterator i=nb_list.begin();
 	 i!=nb_list.end();
 	 i++) {
 	//    cerr << "(" << i->c->index << " " << i->Dij << ")";
+}
 //  cerr << "}" << endl;
+}
 */
 // Tests whether Cell p (given as Vector, remember that Cell is a
 // Vector) is within polygon formed by nearest neighbor cells
 //
 // Based on algorithm and code by Paul Bourke, see
 // http://astronomy.swin.edu.au/~pbourke/geometry/insidepoly/
 //
 // Note: works for 2D only; projects everything on z=0;
 /*
 #define MIN(x,y) (x < y ? x : y)
 #define MAX(x,y) (x > y ? x : y)
  */
 /*bool Cell::CellInsidePolygonP(Vector &p)
+{
     int counter = 0;
     double xinters;
     Vector p1,p2;
     //p1 = polygon[0];
     p1 = *(nb_list.begin()->c);
     int N=nb_list.size();
     list<Neighbor>::const_iterator nb=nb_list.begin();
     for (int i=1;i<=N;i++) {
 		nb++;
 		if (nb!=nb_list.end()) {
 			p2 = *(nb->c);
 		} else {
 			p2 = *(nb_list.begin()->c);
+		}
 		if (p.y > MIN(p1.y,p2.y)) {
 			if (p.y <= MAX(p1.y,p2.y)) {
 				if (p.x <= MAX(p1.x,p2.x)) {
 					if (p1.y != p2.y) {
 						xinters = (p.y-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x;
 						if (p1.x == p2.x || p.x <= xinters)
 							counter++;
+					}
+				}
+			}
+		}
 		p1 = p2;
+    }
     if (counter % 2 == 0)
 		return false;
     else
 		return true;
 }*/
 /* // at new position cell should be able to "see" all polygon sides
 bool Cell::NewPointValidP(Vector &p) {
 	//int ninvtri=0;
 	for (list<Neighbor>::const_iterator nb=nb_list.begin();
 		 nb!=nb_list.end();
 		 nb++) {
 		Vector p1=*(nb->c); // first neighbor
 		list<Neighbor>::const_iterator nextv=nb; nextv++;
 		if (nextv==nb_list.end()) {
 			if (Boundary()==None) {
 				nextv=nb_list.begin();
 			} else continue;
+		}
 		Vector p2=*(nextv->c);
 		Vector v1=(p1-p);
 		Vector v2=(p2-p1);
 		Vector cross=v1*v2;
 		//    //cerr << "[" << cross << "]" << endl;
 		if (cross.z<0) {
 			// One of the triangles has "inverted".
 			//if (Boundary()==None || ninvtri)
 			return false;
 			//else
 			// accept one "inverted" triangle
 			//ninvtri++;
+		}
+	}
 	return true;
 }*/
 // void Cell::CheckForDivision(void) {
 //  //  if (/* Chemical(0)<0.4 && */ /* differentiated cells do not divide */ area > 2*base_area /* || Length()>50 */) {
 //  if (area > par.rel_cell_div_threshold * base_area ) {
 //    /* remark no longer valid? //m->IncreaseCellCapacityIfNecessary();
 //      // Note that calling Divide as follows prevents trouble if cell
 //      // vector is relocated */
 //      Divide();
 //  }
 //}
 /* void Cell::CheckForGFDrivenDivision(void) {
 if (area > base_area && chem[0]>par.gf_div_threshold) {
     //int ind=index;
     if (index==1) return; // petiole does not divide
     // remark no longer valid?
     //m->IncreaseCellCapacityIfNecessary();
     // Note that calling Divide as follows prevents trouble if cell
     // vector is relocated
     Vector horizontal(1,0,0);
     Vector vertical(0,1,0);
     double r;
     if ((r=RANDOM())>par.vertdivprob) {
 		DivideOverAxis(horizontal);
     } else {
 		cerr << "[" << r << "]";
 		DivideOverAxis(vertical);
+    }
+}
+}
 */
 // return (a measure of) the strain of this cell
 /*Vector CellBase::Strain(void) const {
 	cerr << "Sorry, CellBase::strain currently not implemented" << endl;
 	std::exit(1);
 	// Reason: we do not want to include "Node" in the plugins (Node::target_length below), and we do need Strain anyway...
 	// go over all wall elements of the cell
 	 Vector Fvec;
 	for (list<Node *>::const_iterator n=nodes.begin();
 		 n!=nodes.end();
 		 n++) {
 		list<Node *>::const_iterator nn=n; nn++;
 		if (nn==nodes.end()) nn=nodes.begin();
 		Vector wall_element = *(*n) - *(*nn);
 		// assume k=1 (Hooke's constant), for now
 		double Fscal  = (Node::target_length - wall_element.Norm())/Node::target_length;
 		Fvec += Fscal * wall_element.Normalised();
+	}
 	return Fvec;
 } */
 /* void Cell::Flux(double *flux, double D)  {
 // Algorithm according to Rudge & Haseloff 2005
 // (will we need to take cell area into account?)
 // For the time being, we don't: assume cell area is
 // mainly determined by vacuole.
 // Currently implements Rolland-Lagan-Mitchison algorithm
 // Rolland-Lagan and Prusinkiewicz, The Plant Journal (2005), 44, 854-865
 // currently I only implemented passive, diffusive transport
 // active transport will be added later
 // loop over cell edges
 for (int c=0;c<Cell::NChem();c++) flux[c]=0.;
 for (list<Wall>::iterator i=walls.begin();
 	 i!=walls.end();
 	 i++) {
 	// leaf cannot take up chemicals from environment ("no flux boundary")
 	if (i->c2 < 0) continue;
 	// calculate edge length
 	// (will later be updated during node displacement for efficiency)
 	double edge_length = (m->nodes[i->n1]-m->nodes[i->n2]).Norm();
 	// D is "background diffusion coefficient" (Rolland-Lagan)
 	// flux depends on edge length and concentration difference */
 	// i->phi = edge_length * ( /* i->D +*/ D ) * ( m->cells[i->c2].chem[0] - chem[0] );
 	/*
 	 if (m->cells[i->c1].index!=index) {
 		 cerr << "Warning, bad cells boundary: " << m->cells[i->c1].index << ", " << index << endl;
+	 }
 	 flux[0] += i->phi;
 	 //double deltaD = par.alpha * (i->phi*i->phi) - par.gamma * i->D; // Note beta=0
 	 //i->D += par.dt*deltaD;
 	 //cerr << "[ i->D = " << i->D << ", deltaD = " << deltaD << "]";
 	 //if (i->D > par.Dmax) i->D=par.Dmax;
 	 // first calculate all fluxes, we update diffusion coefficient afterwards.
 	 // cerr << "[ " << edge_length << ", " << m->cells[i->c2].chem[0] << " - " << chem[0] << "]";
+}
+}
 */
 	// Save the cell to a stream so we can reconstruct its state later
 	void CellBase::Dump(ostream &os) const {
 		os << index << " " << nodes.size() << endl;
 		Vector::Dump(os);
 		os << endl;
 		for (list<Node *>::const_iterator i=nodes.begin();i!=nodes.end();i++) {
 			os << *i << " ";
+		}
 		os << endl;
 		os << index << " " << neighbors.size() << endl;
 		for (list<CellBase *>::const_iterator i=neighbors.begin();i!=neighbors.end();i++) {
 			os << *i << " ";
+		}
 		os << endl;
 		os << walls.size() << endl;
 		/*for (list<Wall *>::const_iterator i=walls.begin();i!=walls.end(); i++) {
 			(*i)->Dump(os);
 		}*/
 		os << endl;
 		os << NChem() << " ";
 		for (int i=0;i<NChem();i++) {
 			os << chem[i] << " ";
+		}
 		os << endl;
 		os << NChem() << " ";
 		for (int i=0;i<NChem();i++) {
 			os << new_chem[i] << " ";
+		}
 		os << endl;
 		os << boundary << " " << area << " " << target_area << " " << target_length << " " << fixed << " " << intgrl_xx << " " << intgrl_xy << " " << intgrl_yy << " " << intgrl_x << " " << intgrl_y << " " << source << " ";
 		cellvec.Dump(os);
 		os << " " << source_conc << " " << source_chem;
 		os << endl;
+	}
 	void CellBase::UnfixNodes(void) {
 		for (list<Node *>::const_iterator i=nodes.begin();
 			 i!=nodes.end();
 			 i++) {
 			(*i)->Unfix();
+		}
+	}

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;
+	}
     //Vector Strain(void) const;
 	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 );
 	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);
 inline Vector PINdir(CellBase *here, CellBase *nb, Wall *w) {
 	return w->getTransporter( here, 1)  *  w->getInfluxVector(here);
+}
 #endif

src/libplugin.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 += staticlib
 QMAKE_CXXFLAGS += -fexceptions
 QMAKE_CXXFLAGS_DEBUG += -g3
 QMAKE_CXXFLAGS_DEBUG += -DQDEBUG
 QMAKE_CXXFLAGS_DEBUG -= -finstrument-functions
 DEFINES = QTGRAPHICS # VLEAFPLUGIN
-DESTDIR = build_models/lib
+DESTDIR = ../lib
 PERLDIR = ./perl
 PARTMPL = VirtualLeafpar.tmpl
 QT += qt3support
 TARGET = vleaf
 TEMPLATE = lib
 HEADERS = \
  cellbase.h \
  matrix.h \
  output.h \
  parameter.h \
  parse.h \
  random.h \
  simplugin.h \
  UniqueMessage.h \
  vector.h \
  wallbase.h \
  warning.h
 SOURCES = \
  cellbase.cpp \
  matrix.cpp \
  output.cpp \
  parameter.cpp \
  parse.cpp \
  random.cpp \
  simplugin.cpp \
  UniqueMessage.cpp \
  vector.cpp \
  wallbase.cpp \
  warning.cpp
 unix {
  system(rm -f parameter.cpp parameter.h) # this is performed here when qmake is envoked and not in the resulting makefile.
  QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2
  QMAKE_LFLAGS += -fPIC
  LIBS += -lxml2 -lz -lm
+}
 win32 {
  system(DEL parameter.cpp parameter.h)
  LIBXML2DIR = C:\libxml2
  LIBICONVDIR = C:\libiconv
  LIBZDIR = C:\libz
  QMAKE_CXXFLAGS += -DLIBXML_STATIC
  QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include
+}
 system(perl $$PERLDIR/make_parameter_source.pl $$PARTMPL)
 system(perl $$PERLDIR/make_pardialog_source.pl $$PARTMPL)
 #system(perl $$PERLDIR/make_xmlwritecode.pl -h $$REACTIONS)
+#

src/mesh.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 <string>
 #include <algorithm>
 #include <vector>
 #include <sstream>
 #include <cstdlib>
 //#include <cerrno>
 #include <cstring>
 #include <numeric>
 #include <functional>
 #include <fstream>
 #include <QPair>
 #include "mesh.h"
 #include "tiny.h"
 #include "parameter.h"
 #include "random.h"
 #include "pi.h"
 #include "parse.h"
 #include "matrix.h"
 #include "sqr.h"
 #include "nodeset.h"
 #include "nodeitem.h"
 #include "simplugin.h"
 #include <QDebug>
 #include <set>
 #include <iostream>
 #include <iterator>
 static const std::string _module_id("$Id$");
+static const std::string _module_id("$Id: mesh.cpp,v 79f94eaa3b9e 2010/04/14 07:51:02 michael $");
 extern Parameter par;
 void Mesh::AddNodeToCellAtIndex(Cell *c, Node *n, Node *nb1, Node *nb2, list<Node *>::iterator ins_pos) {
   c->nodes.insert(ins_pos, n);
   n->owners.push_back( Neighbor(c, nb1, nb2 ) );
+}
 void Mesh::AddNodeToCell(Cell *c, Node *n, Node *nb1, Node *nb2) {
   c->nodes.push_back( n );
   n->owners.push_back( Neighbor(c, nb1, nb2 ) );
+}
 void Mesh::PerturbChem(int chemnum, double range) {
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
     (*i)->chem[chemnum] += range*(RANDOM()-0.5);
     if ((*i)->chem[chemnum]<0.) (*i)->chem[chemnum]=0.;
+  }
+}
 void Mesh::CellFiles(const Vector ll, const Vector ur) {
   Cell *cell = RectangularCell(ll,ur,0.001);
   for (int c=0;c<Cell::NChem();c++) {
     cell->SetChemical(c,par.initval[c]);
+  }
   cell->SetTargetArea(cell->CalcArea());
   Vector axis(1,0,0);
   // divide rectangle a number of times
   for (int i=0;i<6;i++) {
     IncreaseCellCapacityIfNecessary();
     vector <Cell *> current_cells = cells;
     for (vector<Cell *>::iterator j=current_cells.begin();
 	 j!=current_cells.end();j++) {
       (*j)->DivideOverAxis(axis);
+    }
     axis=axis.Perp2D();
+  }
   IncreaseCellCapacityIfNecessary();
   axis=axis.Perp2D();
   vector <Cell *> current_cells = cells;
   for (vector<Cell *>::iterator j=current_cells.begin();
 	 j!=current_cells.end();j++) {
    (*j)->DivideOverAxis(axis);
+ }
  double sum_l=0; int n_l=0;
   for (list<Node *>::const_iterator i=cell->nodes.begin();
        i!=cell->nodes.end();
        i++) {
     list<Node *>::const_iterator nb=i; nb++;
     if (nb==cell->nodes.end())
       nb=cell->nodes.begin();
     double l = (**nb-**i).Norm();
     sum_l += l;
     n_l++;
+  }
   Node::target_length = sum_l/(double)n_l;
   // a bit more tension
   Node::target_length/=4.;
   SetBaseArea();
+}
 Cell *Mesh::RectangularCell(const Vector ll, const Vector ur, double rotation) {
   Cell *cell=AddCell(new Cell());
   cell->m=this;
   Matrix rotmat;
   rotmat.Rot2D(rotation); // rotation over 0,0
   Node *n1=AddNode(new Node(rotmat * ll));
   Node *n2=AddNode(new Node(rotmat * Vector(ll.x, ur.y,0)));
   Node *n3=AddNode(new Node(rotmat * ur));
   Node *n4=AddNode(new Node(rotmat * Vector(ur.x, ll.y,0)));
   n1->boundary=true;
   n2->boundary=true;
   n3->boundary=true;
   n4->boundary=true;
   //n1.fixed=n2.fixed=n3.fixed=n4.fixed=true;
   AddNodeToCell(cell, n4,
 		n1,
 		n3);
   AddNodeToCell(cell, n3,
 		n4,
 		n2);
   AddNodeToCell(cell, n2,
 		n3,
 		n1);
   AddNodeToCell(cell, n1,
 		n2,
 		n4);
   AddNodeToCell(boundary_polygon, n4,
 		n1,
 		n3);
   AddNodeToCell(boundary_polygon, n3,
 		n4,
 		n2);
   AddNodeToCell(boundary_polygon, n2,
 		n3,
 		n1);
   AddNodeToCell(boundary_polygon, n1,
 		n2,
 		n4);
   cell->setCellVec(Vector(0,1,0));
   boundary_polygon->m = this;
   boundary_polygon->area = 0;
   cell->area = cell->CalcArea();
   // target length is the length of the elements
   Node::target_length = ur.y-ll.y;
   // a bit more tension
   Node::target_length/=2;
   cell->SetIntegrals();
   cell->ConstructNeighborList();
   return cell;
+}
 Cell &Mesh::EllipticCell(double xc, double yc, double ra, double rb,  int nnodes, double rotation) {
   int first_node=Node::nnodes;
   //  nodes.reserve(nodes.size()+nnodes);
   //cells.push_back(Cell(xc,yc));
   Cell *c=AddCell(new Cell(xc,yc));
   c->m=this;
   for (int i=0;i<nnodes;i++) {
     double angle=2*Pi*(i/(double)nnodes);
     double x=xc+ra*cos(angle)*cos(rotation) - rb*sin(angle)*sin(rotation);
     double y=yc+ra*cos(angle)*sin(rotation) + rb*sin(angle)*cos(rotation);
     Node *n=AddNode(new Node(x,y,0));
     n->boundary = true;
+  }
   for (int i=0;i<nnodes;i++) {
     AddNodeToCell(c,
 		  nodes[first_node + i],
 		  nodes[first_node+ (nnodes+i-1)%nnodes],
 		  nodes[first_node+ (i + 1)%nnodes]);
     AddNodeToCell(boundary_polygon,
 		  nodes[first_node + i],
 		  nodes[first_node+ (nnodes+i-1)%nnodes],
 		  nodes[first_node+ (i + 1)%nnodes]);
+  }
   boundary_polygon->m = this;
   boundary_polygon->area = 0;
   c->area = c->CalcArea();
   // target length is the length of the elements
   Node::target_length = (2 * ((ra +rb)/2.) * sin (Pi/nnodes));
   // a bit more tension
   Node::target_length/=2;
   //boundary_polygon = c;
   /*  list<int>::iterator nb;
   for (list<int>::iterator i=c->nodes.begin();
        i!=c->nodes.end();
        i++) {
     nb = i; nb++;
     if (nb==c->nodes.end()) {
       nb=c->nodes.begin();
+    }
     int next = *nb;
     nb = i;
     if (nb==c->nodes.begin()) {
       nb=c->nodes.end();
+    }
     nb--;
     int previous = *nb;
     getNode(*i).cells.push_back( Neighbor(boundary_polygon->index, next, previous) );
   }*/
   c->SetIntegrals();
   //c->ConstructNeighborList();
   //c->ConstructWalls();
   // initial cell has one wall with the outside world
   //c->walls.push_back( new Wall ( nodes.front(), nodes.front(), c, boundary_polygon ));
   c->at_boundary=true;
   return *c;
+}
 Cell &Mesh::LeafPrimordium(int nnodes, double pet_length) {
   // first leaf cell
   int first_node=Node::nnodes;
   Cell *circle=AddCell(new Cell(0,0));
   circle->m=this;
   const double ra=10, rb=10;
   const double xc=0,yc=0;
   const double rotation=0;
   for (int i=0;i<nnodes;i++) {
     double angle=2*Pi*(i/(double)nnodes);
     double x=xc+ra*cos(angle)*cos(rotation) - rb*sin(angle)*sin(rotation);
     double y=yc+ra*cos(angle)*sin(rotation) + rb*sin(angle)*cos(rotation);
     Node *n=AddNode(new Node(x,y,0));
     /* if (angle > 1.25*Pi && angle < 1.75*Pi ) {
       n.sam = true;
       }*/
     AddNodeToCell(circle,
 		  n,
 		  nodes[first_node+ (nnodes+i-1)%nnodes],
 		  nodes[first_node+ (i + 1)%nnodes]);
+  }
   boundary_polygon->m = this;
   boundary_polygon->area = 0;
   circle->area = circle->CalcArea();
   // target length is the length of the elements
   Node::target_length = (2 * ((ra +rb)/2.) * sin (Pi/nnodes));
   // a bit more tension
   Node::target_length/=2;
   circle->SetIntegrals();
   //return c;
   circle->SetTargetArea(2*circle->Area());
   // Petiole: starts at both sides of the circular cell
   // get position of the (n/4)'th and (3*(n/4))'th node.
   list<Node *>::reverse_iterator it_n1=circle->nodes.rbegin();
   for (int i=0;i<nnodes/2;i++)
     it_n1++;
   it_n1--;
   list<Node *>::reverse_iterator it_n2=--circle->nodes.rend();
   /* for (int i=0;i<n/2;i++)
      it_n2++;*/
   Cell *petiole=AddCell(new Cell());
   Node *n1 = *it_n1;
   Node *n2 = *it_n2;
   Node *n3=AddNode( new Node ( *n2 + Vector( 0, pet_length, 0) ) );
   Node *n4=AddNode( new Node ( *n1 + Vector( 0, pet_length, 0) ) );
   n3->boundary=true;
   n4->boundary=true;
   AddNodeToCell(petiole, *it_n1,
 		n4,
 		nodes[(*it_n2)->Index()
 		      + (( (*it_n1)->Index() - (*it_n2)->Index() )-1+nnodes)%nnodes]);
   list<Node *>::reverse_iterator i=it_n1; i++;
   for (;
        i!=it_n2;
        //(++i) == circle->nodes.rend() ? i : i=circle->nodes.rbegin() ) {
        ++i) {
     AddNodeToCell(petiole,
 		  *i,
 		  nodes[(*it_n2)->Index() + (((*i)->Index()-(*it_n2)->Index()) + 1)%nnodes],
 		  nodes[(*it_n2)->Index() + (((*i)->Index()-(*it_n2)->Index())-1+nnodes)%nnodes]);
+  }
   AddNodeToCell(petiole, *it_n2,
 		*it_n2 + 1,
 		n3);
   (*it_n2)->boundary=true;
   //petiole.nodes.push_back(n3.Index());
   //petiole.nodes.push_back(n4.Index());
   AddNodeToCell(petiole,
 		n3,
 		n2,
 		n4);
   AddNodeToCell(petiole,
 		n4,
 		n3,
 		n1);
   cerr << circle << endl;
   cerr << petiole << endl;
   AddNodeToCell(boundary_polygon, *it_n1,
 		n4,
 		*it_n2 + ((*it_n1-*it_n2)+1)%nnodes); // is this gonna work?
   (*it_n1)->boundary=true;
   for (int i=0;i<nnodes;i++) {
     if (nodes[(first_node + i)]->owners.size()==1) {
       AddNodeToCell(boundary_polygon,
 		    nodes[first_node +i],
 		    nodes[first_node+ (nnodes+i-1)%nnodes],
 		    nodes[first_node+ (i + 1)%nnodes]);
       nodes[first_node+i]->boundary=true;
+    }
+  }
   AddNodeToCell(boundary_polygon, *it_n2,
 		nodes[(nnodes+(*it_n2)->Index() - 1)%nnodes],
 		n3);
   AddNodeToCell(boundary_polygon,
 		n3,
 		n2,
 		n4);
   AddNodeToCell(boundary_polygon,
 		n4,
 		n3,
 		n1);
   // make petiole solid
   for (list<Node *>::iterator i=petiole->nodes.begin();
        i!=petiole->nodes.end();
        i++) {
     (*i)->Fix();
+  }
   petiole->Fix();
   petiole->area=petiole->CalcArea();
   petiole->target_area=petiole->area;
   petiole->ConstructNeighborList();
   circle->ConstructNeighborList();
   boundary_polygon->ConstructConnections();
   boundary_polygon->ConstructNeighborList();
   circle->setCellVec(Vector(0,1,0));
   return *circle;
+}
 /*Cell &Mesh::Box() {
 }*/
 // return bounding box of mesh
 void Mesh::BoundingBox(Vector &LowerLeft, Vector &UpperRight) {
   LowerLeft = **nodes.begin();
   UpperRight = **nodes.begin();
   for (vector<Node *>::iterator c=nodes.begin();
        c!=nodes.end();
        c++) {
     if ((*c)->x < LowerLeft.x)
       LowerLeft.x = (*c)->x;
     if ((*c)->y < LowerLeft.y)
       LowerLeft.y = (*c)->y;
     if ((*c)->z < LowerLeft.z)
       LowerLeft.z = (*c)->z;
     if ((*c)->x > UpperRight.x)
       UpperRight.x = (*c)->x;
     if ((*c)->y > UpperRight.y)
       UpperRight.y = (*c)->y;
     if ((*c)->z > UpperRight.z)
       UpperRight.z = (*c)->z;
+  }
+}
 double Mesh::Area(void) {
   double area=0;
   vector<Cell *>::iterator i=cells.begin();
   while (i != cells.end()) {
     area += (*(i++))->Area();
+  }
   return area;
+}
 void Mesh::SetBaseArea(void) {
   // Set base area to mean area.
   // This method is typically called during initiation, after
   // defining the first cell
 	Cell::BaseArea() = Area()/cells.size();
+}
 // for optimization, we moved Displace to Mesh
 class DeltaIntgrl {
 public:
   double area;
   double ix, iy;
   double ixx,ixy,iyy;
   DeltaIntgrl(double sarea,double six,double siy,double sixx,double sixy,double siyy) {
     area=sarea;
     ix=six;
     iy=siy;
     ixx=sixx;
     ixy=sixy;
     iyy=siyy;
+  }
 };
 void Mesh::Clear(void) {
 	// clear nodes
 	for (vector<Node *>::iterator i=nodes.begin();
 		 i!=nodes.end();
 		 i++) {
 		delete *i;
+	}
 	nodes.clear();
 	Node::nnodes=0;
 	node_insertion_queue.clear();
 	// Clear NodeSets
 	for (vector<NodeSet *>::iterator i=node_sets.begin();
 		 i!=node_sets.end();
 		 i++) {
 		delete *i;
+	}
 	node_sets.clear();
 	time = 0;
 	// clear cells
 	for (vector<Cell *>::iterator i=cells.begin();
 		 i!=cells.end();
 		 i++) {
 		delete *i;
+	}
 	cells.clear();
 	Cell::NCells() = 0;
 	delete boundary_polygon;
 	// Clear walls
 	for (list<Wall *>::iterator i=walls.begin();
 		 i!=walls.end();
 		 i++) {
 		delete *i;
+	}
 	walls.clear();
 	WallBase::nwalls = 0;
 	//tmp_walls->clear();
 	shuffled_cells.clear();
 	shuffled_nodes.clear();
 	//Cell::ncells=0;
 	/*    Cell::ClearNCells();
 	 Node::nnodes=0;
 	 cells.clear();
 	 nodes.clear();
 	 shuffled_cells.clear();
 	 shuffled_nodes.clear();
 	 node_insertion_queue.empty();
 	 cerr << "Meshed cleared: cells: " << cells.size() << ", nodes: " << nodes.size() << endl;
 	 */
 	cerr << "cells.size() = " << cells.size() << endl;
 	cerr << "walls.size() = " << walls.size() << endl;
 	cerr << "nodes.size() = " << nodes.size() << endl;
+}
 double Mesh::DisplaceNodes(void) {
   MyUrand r(shuffled_nodes.size());
   random_shuffle(shuffled_nodes.begin(),shuffled_nodes.end(),r);
   double sum_dh=0;
   list<DeltaIntgrl> delta_intgrl_list;
   for_each( node_sets.begin(), node_sets.end(), mem_fun( &NodeSet::ResetDone ) );
   for (vector<Node *>::const_iterator i=shuffled_nodes.begin();
        i!=shuffled_nodes.end();
        i++) {
     //int n=shuffled_nodes[*i];
     Node &node(**i);
     // Do not allow displacement if fixed
     //if (node.fixed) continue;
     if (node.DeadP()) continue;
     // Attempt to move this cell in a random direction
     double rx=par.mc_stepsize*(RANDOM()-0.5); // was 100.
     double ry=par.mc_stepsize*(RANDOM()-0.5);
     // Uniform with a circle of radius par.mc_stepsize
     /* double r = RANDOM() * par.mc_stepsize;
        double th = RANDOM()*2*Pi;
        double rx = r * cos(th);
        double ry = r * sin(th);
     */
     Vector new_p(node.x+rx,node.y+ry,0);
     Vector old_p(node.x,node.y,0);
     /* if (node.boundary  && boundary_polygon->MoveSelfIntersectsP(n,  new_p )) {
     // reject if move of boundary results in self intersection
     continue;
     }*/
     if (node.node_set) {
 		// move each node set only once
 		if (!node.node_set->DoneP())
 			node.node_set->AttemptMove(rx,ry);
     } else {
       // for all cells to which this node belongs:
       //   calculate energy difference
       double area_dh=0.;
       double length_dh=0.;
       double bending_dh=0.;
       double cell_length_dh=0.;
       double alignment_dh=0.;
       double old_l1=0.,old_l2=0.,new_l1=0.,new_l2=0.;
       double sum_stiff=0.;
       double dh=0.;
       for (list<Neighbor>::const_iterator cit=node.owners.begin();
 	   cit!=node.owners.end();
 	   cit++) {
 	//Cell &c=m->getCell(cit->cell);
 	//Cell &c=cit->cell->BoundaryPolP()?*boundary_polygon:*(cit->cell);
 	Cell &c=*((Cell *)(cit->cell));
 	//Cell &c=cells[cit->cell];
 	if (c.MoveSelfIntersectsP(&node,  new_p )) {
 	  // reject if move results in self intersection
 	  //
 	  // I know: using goto's is bad practice... except when jumping out
 	  // of deeply nested loops :-)
 	  //cerr << "Rejecting due to self-intersection\n";
 	  goto next_node;
+	}
 	// summing stiffnesses of cells. Move has to overcome this minimum required energy.
 	sum_stiff += c.stiffness;
 	// area - (area after displacement): see notes for derivation
 	//Vector i_min_1 = m->getNode(cit->nb1);
 	Vector i_min_1 = *(cit->nb1);
       //Vector i_plus_1 = m->getNode(cit->nb2);
 	Vector i_plus_1 = *(cit->nb2);
 	// We must double the weights for the perimeter (otherwise they start bulging...)
 	double w1, w2;
 	if (node.boundary && cit->nb1->boundary)
 #ifdef FLEMING
 	  w1 = par.rel_perimeter_stiffness;
 #else
 	w1=2;
 #endif
 	else
 	  w1 = 1;
 	if (node.boundary && cit->nb2->boundary)
 #ifdef FLEMING
 	  w2 = par.rel_perimeter_stiffness;
 #else
 	w2 = 2;
 #endif
 	else
 	  w2 = 1;
 	//if (cit->cell>=0) {
 	if (!cit->cell->BoundaryPolP()) {
 	  double delta_A = 0.5 * ( ( new_p.x - old_p.x ) * (i_min_1.y - i_plus_1.y) +
 				   ( new_p.y - old_p.y ) * ( i_plus_1.x - i_min_1.x ) );
 	  area_dh +=  delta_A * (2 * c.target_area - 2 * c.area + delta_A);
 	  // cell length constraint
 	  // expensive and not always needed
 	  // so we check the value of lambda_celllength
 	  if (/* par.lambda_celllength */  cit->cell->lambda_celllength) {
 	    double delta_ix =
 	      (i_min_1.x + new_p.x)
 	      * (new_p.x * i_min_1.y- i_min_1.x * new_p.y) +
 	      (new_p.x + i_plus_1.x)
 	      * (i_plus_1.x * new_p.y- new_p.x * i_plus_1.y) -
 	      (i_min_1.x + old_p.x)
 	      * (old_p.x * i_min_1.y- i_min_1.x * old_p.y) -
 	      (old_p.x + i_plus_1.x)
 	      * (i_plus_1.x * old_p.y - old_p.x * i_plus_1.y);
 	    double delta_iy =
 	      (i_min_1.y + new_p.y)
 	      * (new_p.x * i_min_1.y- i_min_1.x * new_p.y) +
 	      (new_p.y + i_plus_1.y)
 	      * (i_plus_1.x * new_p.y- new_p.x * i_plus_1.y) -
 	      (i_min_1.y + old_p.y)
 	      * (old_p.x * i_min_1.y- i_min_1.x * old_p.y) -
 	      (old_p.y + i_plus_1.y)
 	      * (i_plus_1.x * old_p.y - old_p.x * i_plus_1.y);
 	    double delta_ixx =
 	      (new_p.x*new_p.x+
 	       i_min_1.x*new_p.x+
 	       i_min_1.x*i_min_1.x ) *
 	      (new_p.x*i_min_1.y - i_min_1.x*new_p.y) +
 	      (i_plus_1.x*i_plus_1.x+
 	       new_p.x*i_plus_1.x+
 	       new_p.x*new_p.x ) *
 	      (i_plus_1.x*new_p.y - new_p.x*i_plus_1.y) -
 	      (old_p.x*old_p.x+
 	       i_min_1.x*old_p.x+
 	       i_min_1.x*i_min_1.x ) *
 	      (old_p.x*i_min_1.y - i_min_1.x*old_p.y) -
 	      (i_plus_1.x*i_plus_1.x+
 	       old_p.x*i_plus_1.x+
 	       old_p.x*old_p.x ) *
 	      (i_plus_1.x*old_p.y - old_p.x*i_plus_1.y);
 	    double delta_ixy =
 	      (i_min_1.x*new_p.y-
 	       new_p.x*i_min_1.y)*
 	      (new_p.x*(2*new_p.y+i_min_1.y)+
 	       i_min_1.x*(new_p.y+2*i_min_1.y)) +
 	      (new_p.x*i_plus_1.y-
 	       i_plus_1.x*new_p.y)*
 	      (i_plus_1.x*(2*i_plus_1.y+new_p.y)+
 	       new_p.x*(i_plus_1.y+2*new_p.y)) -
 	      (i_min_1.x*old_p.y-
 	       old_p.x*i_min_1.y)*
 	      (old_p.x*(2*old_p.y+i_min_1.y)+
 	       i_min_1.x*(old_p.y+2*i_min_1.y)) -
 	      (old_p.x*i_plus_1.y-
 	       i_plus_1.x*old_p.y)*
 	      (i_plus_1.x*(2*i_plus_1.y+old_p.y)+
 	       old_p.x*(i_plus_1.y+2*old_p.y));
 	    double delta_iyy =
 	      (new_p.x*i_min_1.y-
 	       i_min_1.x*new_p.y)*
 	      (new_p.y*new_p.y+
 	       i_min_1.y*new_p.y+
 	       i_min_1.y*i_min_1.y ) +
 	      (i_plus_1.x*new_p.y-
 	       new_p.x*i_plus_1.y)*
 	      (i_plus_1.y*i_plus_1.y+
 	       new_p.y*i_plus_1.y+
 	       new_p.y*new_p.y ) -
 	      (old_p.x*i_min_1.y-
 	       i_min_1.x*old_p.y)*
 	      (old_p.y*old_p.y+
 	       i_min_1.y*old_p.y+
 	       i_min_1.y*i_min_1.y ) -
 	      (i_plus_1.x*old_p.y-
 	       old_p.x*i_plus_1.y)*
 	      (i_plus_1.y*i_plus_1.y+
 	       old_p.y*i_plus_1.y+
 	       old_p.y*old_p.y );
 	    delta_intgrl_list.push_back(DeltaIntgrl(delta_A,delta_ix,delta_iy,delta_ixx,delta_ixy,delta_iyy));
 	    Vector old_axis;
 	    double old_celllength = c.Length(&old_axis);
 	    old_axis=old_axis.Normalised().Perp2D();
 	    // calculate length after proposed update
 	    double intrx=(c.intgrl_x-delta_ix)/6.;
 	    double intry=(c.intgrl_y-delta_iy)/6.;
 	    double ixx=((c.intgrl_xx-delta_ixx)/12.)-(intrx*intrx)/(c.area-delta_A);
 	    double ixy=((c.intgrl_xy-delta_ixy)/24.)+(intrx*intry)/(c.area-delta_A);
 	    double iyy=((c.intgrl_yy-delta_iyy)/12.)-(intry*intry)/(c.area-delta_A);
 	    double rhs1=(ixx+iyy)/2., rhs2=sqrt( (ixx-iyy)*(ixx-iyy)+4*ixy*ixy )/2.;
@@ @@ -1218,1049 +1218,1061 @@ void Mesh::CircumCircle(double x1,double @@
     *yc = m1 * (*xc - mx1) + my1;
   } else {
     m1 = - (x2-x1) / (y2-y1);
     m2 = - (x3-x2) / (y3-y2);
     mx1 = (x1 + x2) / 2.0;
     mx2 = (x2 + x3) / 2.0;
     my1 = (y1 + y2) / 2.0;
     my2 = (y2 + y3) / 2.0;
     *xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
     *yc = m1 * (*xc - mx1) + my1;
+  }
   dx = x2 - *xc;
   dy = y2 - *yc;
   rsqr = dx*dx + dy*dy;
   *r = sqrt(rsqr);
   return;
   // Suggested
   // return((drsqr <= rsqr + EPSILON) ? TRUE : FALSE);
+}
 //
 // return the total amount of chemical "ch" in the leaf
 double Mesh::SumChemical(int ch) {
   double sum=0.;
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
     sum+=(*i)->chem[ch];
+  }
   return sum;
+}
 void Mesh::CleanUpCellNodeLists(void) {
   typedef vector <vector<Cell *>::iterator> CellItVect;
   CellItVect cellstoberemoved;
   vector<int> cellind;
   // Start of by removing all stale walls.
   //DeleteLooseWalls();
   // collect all dead cells that need to be removed from the simulation
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
     if ((*i)->DeadP()) {
       // collect the iterators
       cellstoberemoved.push_back(i);
       // collect the indices
       cellind.push_back((*i)->index);
     } else {
       // Remove pointers to dead Walls
       for (list<Wall *>::iterator w=(*i)->walls.begin();
 	   w!=(*i)->walls.end();
 	   w++) {
 	if ((*w)->DeadP()) {
 	  (*w)=0;
+	}
+      }
       (*i)->walls.remove(0);
+    }
+  }
   // Remove pointers to dead Walls from BoundaryPolygon
   for (list<Wall *>::iterator w=boundary_polygon->walls.begin();
        w!=boundary_polygon->walls.end();
        w++) {
     if ((*w)->DeadP()) {
       (*w)=0;
+    }
+  }
   boundary_polygon->walls.remove(0);
   // Renumber cells; this is most efficient if the list of dead cell indices is sorted
   sort(cellind.begin(),cellind.end());
   // Reindexing of Cells
   for (vector<int>::reverse_iterator j=cellind.rbegin();
        j!=cellind.rend();
        j++) {
     for (vector<Cell *>::reverse_iterator i=cells.rbegin();
 	 i!=cells.rend();
 	 i++) {
       if (*j < (*i)->index) (*i)->index--;
+    }
+  }
   // Actual deleting of Cells
   // We must delete in reverse order, otherwise the iterators become redefined
   for ( CellItVect::reverse_iterator i=cellstoberemoved.rbegin();
 	i!=cellstoberemoved.rend();
 	i++) {
     Cell::NCells()--;
     cells.erase(*i);
+  }
   // same for nodes
   typedef vector <vector<Node *>::iterator> NodeItVect;
   NodeItVect nodestoberemoved;
   vector<int> nodeindlist;
   // collect iterators and indices of dead nodes
   for (vector<Node *>::iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
     if ((*i)->DeadP()) {
       nodestoberemoved.push_back( i );
       nodeindlist.push_back((*i)->index);
+    }
+  }
   // sort the list of dead nodes for renumbering
   sort(nodeindlist.begin(),nodeindlist.end());
   // Reindicing of Nodes
   for (vector<int>::reverse_iterator j=nodeindlist.rbegin();
        j!=nodeindlist.rend();
        j++) {
     for (vector<Node *>::reverse_iterator i=nodes.rbegin();
 	 i!=nodes.rend();
 	 i++) {
       if (*j < (*i)->index) {
 	(*i)->index--;
+      }
+    }
+  }
   // Actual deleting of nodes
   // We must delete in reverse order, otherwise the iterators become redefined
   for ( NodeItVect::reverse_iterator i=nodestoberemoved.rbegin();
 	i!=nodestoberemoved.rend();
 	i++) {
     Node::nnodes--;
     nodes.erase(*i);
+  }
   for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
     if ((*w)->DeadP()) {
       Wall::nwalls--;
       delete *w;
       *w = 0;
+    }
+  }
   walls.remove( 0 );
   // Clean up all intercellular connections and redo everything
   for (vector<Node *>::iterator i=nodes.begin();
        i!=nodes.end();
        i++) {
     (*i)->owners.clear();
+  }
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
     (*i)->ConstructConnections();
+  }
   boundary_polygon->ConstructConnections();
   /* for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
     delete *w;
+  }
   walls.clear();
   cerr << "Cleared walls\n";
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
     (*i)->ConstructWalls();
+  }
   */
   // remake shuffled_nodes and shuffled cells
   shuffled_nodes.clear();
   shuffled_nodes = nodes;
   shuffled_cells.clear();
   shuffled_cells = cells;
+}
 void Mesh::CutAwayBelowLine( Vector startpoint, Vector endpoint) {
   // Kills all cells below the line startpoint -> endpoint
   Vector perp = (endpoint-startpoint).Perp2D().Normalised();
   cerr << "Before Apoptose\n";
   TestIllegalWalls();
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
     // do some vector geometry to check whether the cell is below the cutting line
     Vector cellvec = ((*i)->Centroid()-startpoint);
     if ( InnerProduct(perp, cellvec) < 0 ) {
       // remove those cells
       (*i)->Apoptose();
+    }
+  }
   cerr << "Before CleanUpCellNodeLists\n";
   TestIllegalWalls();
   CleanUpCellNodeLists();
+}
 void Mesh::CutAwaySAM(void) {
   for (vector<Cell *>::iterator i=cells.begin();
        i!=cells.end();
        i++) {
     if( (*i)->Boundary() == Cell::SAM ) {
 	(*i)->Apoptose();
+    }
+  }
   TestIllegalWalls();
   CleanUpCellNodeLists();
+}
 void Mesh::TestIllegalWalls(void) {
   for (list<Wall *>::iterator w = walls.begin();
        w!=walls.end();
        w++) {
     if ((*w)->IllegalP() ) {
       cerr << "Wall " << **w << " is illegal.\n";
+    }
+  }
+}
 class node_owners_eq : public unary_function<Node, bool> {
   int no;
 public:
   explicit node_owners_eq(int nn) { no=nn; }
   bool operator() (const Node &n) const {
     if (n.CellsSize()==1)
       return true;
     else
       return false;
+  }
 };
 void Mesh::RepairBoundaryPolygon(void) {
   // After serious manipulations (e.g. after cutting part off the
   // leaf) repair the boundary polygon. It assumes the cut line has
   // already been marked "boundary" and the original boundary marks
   // were not removed.
   //
   // So, this function just puts boundary nodes into the boundary
   // polygon in the right order; it cannot detect boundaries from
   // scratch.
   Node *boundary_node=0, *next_boundary_node=0, *internal_node;
   set<int> original_boundary_nodes, repaired_boundary_nodes;
   vector<int> difference; // set difference result
   // Step 0: print out current boundary polygon
 #ifdef QDEBUG
   qDebug() << endl << "Original Boundary Polygon node indices: ";
   foreach (Node* node, boundary_polygon->nodes) {
     qDebug() << node->Index() << " " ;
+  }
   qDebug() << endl << endl;
 #endif
   // Step 1a: Create a set containing the current boundary polygon nodes' Indices.
   foreach (Node* node, boundary_polygon->nodes) {
     original_boundary_nodes.insert(node->Index());
+  }
   // Step 1b: remove all nodes from boundary polygon
   boundary_polygon->nodes.clear();
   // Step 2: Remove all references to the boundary polygon from the Mesh's current list of nodes
   foreach (Node* node, nodes) {
     node->Unmark(); // remove marks, we need them to determine if we have closed the circle
     list<Neighbor>::iterator boundary_ref_pos;
     if ((boundary_ref_pos = find_if (node->owners.begin(), node->owners.end(),
 				     bind2nd(mem_fun_ref(&Neighbor::CellEquals), -1))) != node->owners.end()) {
       // i.e. if one of the node's owners is the boundary polygon
       node->owners.erase(boundary_ref_pos); // remove the reference
+    }
+  }
   // Step 3: Search for the first boundary node.  We reconstruct the
   // boundary polygon by moving along the boundary nodes until we've
   // encircled the polygon. Since manually adding nodes may have
   // turned nodes previously along the boundary into internal nodes,
   // we search through all the node until we find first boundary node
   // and proceed from there. If findNextBoundaryNode() returns a node
   // other than the one passed to it, the original node is the first
   // boundary node.
   foreach (Node* node, nodes) {
     if ((findNextBoundaryNode(node))->index != node->index){
       next_boundary_node = node;
       break;
+    }
+  }
   // We have a problem if we arrive here without having found a boundary node.
   if (!next_boundary_node) throw("Cannot find a boundary node!.");
   // Reconstruct the list of boundary polygon nodes.
   do {
     boundary_node = next_boundary_node;
     boundary_node->Mark();
     boundary_polygon->nodes.push_back(boundary_node);
     next_boundary_node = findNextBoundaryNode(boundary_node);
   } while ( !next_boundary_node->Marked() );
   // Create a set containing the reconstructed boundary polygon nodes' Indices.
   for (list<Node *>::iterator it = boundary_polygon->nodes.begin(); it!=boundary_polygon->nodes.end(); ++it) {
     repaired_boundary_nodes.insert((*it)->Index());
+  }
   // Calculate the difference between the original and repaired sets of boundary nodes
   // yielding the set of nodes that are no longer part of the boundary polygon.
   set_difference(original_boundary_nodes.begin(), original_boundary_nodes.end(),
                  repaired_boundary_nodes.begin(), repaired_boundary_nodes.end(), back_inserter(difference));
   // Tell each node in the difference that it's no longer part of the boundary polygon
   vector<Node *>::iterator internal_node_it;
   foreach (int i, difference){
     internal_node_it = find_if (nodes.begin(), nodes.end(), bind2nd(mem_fun(&Node::IndexEquals), i));
     internal_node = *internal_node_it; // dereference the itterator to get to the node pointer
     if (!internal_node) throw("Found a null Node pointer.");
     internal_node->UnsetBoundary();
+  }
   boundary_polygon->ConstructConnections();
   for (list<Wall *>::iterator w=boundary_polygon->walls.begin();
        w!=boundary_polygon->walls.end();
        w++) {
     if ((*w)->DeadP()) {
       (*w)=0;
+    }
+  }
   boundary_polygon->walls.remove(0);
   boundary_polygon->ConstructNeighborList();
 #ifdef QDEBUG
   cerr << "Repaired Boundary Polygon node indices: ";
   foreach (Node* node, boundary_polygon->nodes){
     qDebug() << node->Index() << " " ;
+  }
   qDebug() << endl ;
   #ifdef _undefined_
   qDebug() << "NODES:" << endl;
   foreach(Node* node, nodes) {
     qDebug() << *node;
+  }
   qDebug() << endl;
   qDebug() << "WALLS:" << endl;
   foreach(Wall* wall, walls) {
     qDebug() << *wall;
+  }
   qDebug() << endl;
   qDebug() << "CELLS:" << endl;
   foreach(Cell* cell, cells) {
     qDebug() << *cell;
+  }
   qDebug() << endl;
   #endif
 #endif
+}
 Node* Mesh::findNextBoundaryNode(Node* boundary_node) {
   bool found_next_boundary_node = false;
   Node *next_boundary_node;
   set<int> boundary_node_owners; // This is a list of the current boundary node's owners' Ids
   vector<int> neighborIds; // A list of the current boundary node's owners' 2nd neighbor Ids
   vector<set<int> *>  nodeOwners; // A vector of set pointers where each set contains the owner Ids of the nodes in the neighborIds list.
   vector<int> intersection; // set intersection result
   // The next boundary node is that which has only one owner in common with the current boundary node
   for (list<Neighbor>::iterator it=boundary_node->owners.begin(); it!=boundary_node->owners.end(); ++it) {
     if (it->cell->Index() != -1) boundary_node_owners.insert(it->cell->Index()); // Save each of the current boundary node's owners' Ids - except the boundary polygon
     set<int> *owners = new set<int>; // create a set to hold a 2nd neighbor's owners' Ids
     nodeOwners.push_back(owners);
     neighborIds.push_back(it->nb2->Index());
     foreach(Neighbor neighbor, it->nb2->owners){
       if (neighbor.cell->Index() != -1) owners->insert(neighbor.cell->Index()); // Save second neighbors' owners' Ids - except the boundary polygon
+    }
+  }
   vector<int>::iterator itt = neighborIds.begin();
   vector<set<int> *>::iterator it = nodeOwners.begin();
   #ifdef QDEBUG
   qDebug() << "Boundary node: " <<  boundary_node->Index() << " is owned by the following cells: ";
   foreach (int i, boundary_node_owners){
     qDebug() << i << "  ";
+  }
   qDebug() << endl;
   #endif
   for (; it < nodeOwners.end(); it++, itt++) {
     intersection.clear();
     set_intersection(boundary_node_owners.begin(), boundary_node_owners.end(), (*it)->begin(), (*it)->end(), back_inserter(intersection));
     #ifdef QDEBUG
     qDebug() << "The intersection of the boundary node(" << boundary_node->Index() << ") owners and its 2nd neighbor(" <<  *itt << ") owners is: ";
     foreach (int i, intersection){
       qDebug() << i << "  ";
+    }
     qDebug() << endl;
     #endif
     if (intersection.size() == 1){
       found_next_boundary_node = true;
       vector<Node *>::iterator next_boundary_node_it = find_if (nodes.begin(), nodes.end(), bind2nd(mem_fun(&Node::IndexEquals), *itt));
       next_boundary_node = *next_boundary_node_it; // defeference the itterator to get to the node pointer
       #ifdef QDEBUG
       qDebug() << "The Current boundary node is: " << boundary_node->Index()
 	       << ". The Next boundary node is: " << *itt << ((next_boundary_node->Marked()) ? " Marked" : " Unmarked") << endl << endl;
       #endif
       break;
+    }
+  }
   #ifdef QDEBUG
   if (!found_next_boundary_node) {
     qDebug() << "OOPS! Didn't find the next boundrary node!" << endl;
+  }
   #endif
   return next_boundary_node;
+}
 void Mesh::CleanUpWalls(void) {
   for (list<Wall *>::iterator w=walls.begin();
        w!=walls.end();
        w++) {
     if ((*w)->DeadP()) {
       delete *w;
       (*w)=0;
+    }
+  }
   walls.remove(0);
+}
 void Mesh::Rotate(double angle, Vector center) {
   // Rotate the mesh over the angle "angle", relative to center point "center".
   Matrix rotmat;
   rotmat.Rot2D(angle);
   for (vector<Node *>::iterator n=nodes.begin();
        n!=nodes.end();
        n++) {
     (*n)->setPos ( rotmat * ( *(*n) - center ) + center );
+  }
+}
 void Mesh::PrintWallList( void ) {
    transform ( walls.begin(), walls.end(), ostream_iterator<Wall>(cerr, "\n"), deref_ptr<Wall> );
+}
 #include <QString>
 //#include "forwardeuler.h"
 #include "rungekutta.h"
 class SolveMesh : public RungeKutta {
 private:
   SolveMesh(void);
 public:
 	SolveMesh(Mesh *m_) {
 		m = m_;
 		kmax=0;
 		kount=0;
 		xp=0; yp=0; dxsav=0;
+	}
 protected:
   virtual void derivs(double x, double *y, double *dydx) {
     // set mesh with new values given by ODESolver
     // (we must do this, because only mesh knows the connections
     // between the variables)
     m->setValues(x,y);
     m->Derivatives(dydx);
     /*static int c=0;
        QString fname("derivs%1.dat");
        ofstream of(fname.arg(++c).ascii());
        for (int i=0;i<m->NEqs();i++) {
        of << x << " " << dxdy[i] << endl;
+       }
        of.close();
     */
     //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++) {
 		if ((*n)->Value()>2 && !(*n)->BoundaryP() ) {
 			QVector<qreal> angles = (*n)->NeighbourAngles();
 			int value_vertex = angles.size();
 			for (QVector<qreal>::ConstIterator i=angles.begin();
 				 i!=angles.end();
 				 i++) {
 				anglesvalues += QPair< qreal, int > (*i, value_vertex);
+			}
+		}
+	}
 	return anglesvalues;
+}
 void Mesh::Clean(void) {
         #ifdef QDEBUG
         qDebug() << "Freeing nodes" << endl;
 	#endif
 	for (vector<Node *>::iterator i=nodes.begin();
 		 i!=nodes.end();
 		 i++) {
 		delete *i;
+	}
 	nodes.clear();
 	Node::nnodes=0;
         #ifdef QDEBUG
         qDebug() << "Freeing node sets" << endl;
         #endif
 	for (vector<NodeSet *>::iterator i=node_sets.begin();
 		 i!=node_sets.end();
 		 i++) {
 		delete *i;
+	}
 	node_sets.clear();
         #ifdef QDEBUG
         qDebug() << "Freeing cells" << endl;
         #endif
 	//CellsStaticDatamembers *old_static_data_mem = Cell::GetStaticDataMemberPointer();
 	for (vector<Cell *>::iterator i=cells.begin();
 		 i!=cells.end();
 		 i++) {
 		delete *i;
+	}
 	//Cell::static_data_members = old_static_data_mem;
 	cells.clear();
 	Cell::NCells()=0;
 	delete boundary_polygon; // (already deleted during cleaning of cells?)
         #ifdef QDEBUG
         qDebug() << "Freeing walls" << endl;
         #endif
 	for (list<Wall *>::iterator i=walls.begin();
 		 i!=walls.end();
 		 i++) {
 		delete *i;
+	}
 	walls.clear();
 	Wall::nwalls=0;
 	node_insertion_queue.clear();
 	shuffled_nodes.clear();
 	shuffled_cells.clear();
 	time = 0.0;
+}
 void Mesh::StandardInit(void) {
 	boundary_polygon = new BoundaryPolygon();
 	Cell &circle=CircularCell(0,0,10,10);
 	circle.SetTargetArea(circle.CalcArea());
 	circle.SetTargetLength(par.target_length);
 	circle.SetLambdaLength(par.lambda_celllength);
 	SetBaseArea();
 	// clean up chemicals
 	for (int c=0; c<Cell::NChem(); c++) {
 		circle.SetChemical(c, 0.);
+	}
+}

src/mesh.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.
+ *
  */
 // Cell derives from Vector, where Vector is simply used as a Vertex
 #ifndef _MESH_H_
 #define _MESH_H_
 #include <vector>
 #include <algorithm>
 #include <queue>
 #include <iterator>
 #include <functional>
 #ifdef QTGRAPHICS
 #include <QGraphicsScene>
 #endif
 #include "cell.h"
 #include "node.h"
 #include "simplugin.h"
 #include <QVector>
 #include <QPair>
 using namespace std;
 // new queue which rejects duplicate elements
 template<class T, class C = deque<T> > class unique_queue : public queue<T,C> {
 public:
 	typedef typename C::value_type value_type;
 	// reimplements push: reject element if it exists already
 	void push(const value_type &x) {
 		if (find (queue<T,C>::c.begin(),queue<T,C>::c.end(),x)==queue<T,C>::c.end()) {
 			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);
 			plugin->CellHouseKeeping(*i);
 			// Call functions of Cell that cannot be called from CellBase, including Division
 			if ((*i)->flag_for_divide) {
 				(*i)->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:
 	vector<Node *> shuffled_nodes;
 	vector<Cell *> shuffled_cells;
 	unique_queue<Edge> node_insertion_queue;
 	BoundaryPolygon *boundary_polygon;
 	double time;
 	SimPluginInterface *plugin;
 	// Private member functions
 	void AddNodeToCell(Cell *c, Node *n, Node *nb1 , Node *nb2);
 	void AddNodeToCellAtIndex(Cell *c, Node *n, Node *nb1 , Node *nb2, list<Node *>::iterator ins_pos);
 	void InsertNode(Edge &e);
 	inline Node *AddNode(Node *n) {
 		nodes.push_back(n);
 		shuffled_nodes.push_back(n);
 		n->m=this;
 		return n;
+	}
 	inline Cell *AddCell(Cell *c) {
 		cells.push_back(c);
 		shuffled_cells.push_back(c);
 		//cerr << "Shuffled cell indices:  ";
 		/*copy(shuffled_cells.begin(),shuffled_cells.end(),ostream_iterator<int>(cerr," "));
 		 cerr << endl;*/
 		c->m=this;
 		return c;
+	}
 	//int Delaunay(void);
 	void CircumCircle(double x1,double y1,double x2,double y2,double x3,double y3,
 					  double *xc,double *yc,double *r);
 	// void RenumberCells(void);
 };
 #endif

src/modelcatalogue.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 <string>
 #include "modelcatalogue.h"
 #include <QVariant>
 static const std::string _module_id("$Id$");
 ModelCatalogue::ModelCatalogue(Mesh *_mesh, Main *_mainwin, const char *model=0) {
 	mesh = _mesh;
 	mainwin = _mainwin;
 	if (model) {
 	  cerr << "Loading model: " << model << endl;
 		LoadPlugin(model);
 	} else {
 	  cerr << "Loading all models." << endl;
 		LoadPlugins();
+	}
+}
 void ModelCatalogue::LoadPlugins() {
 	QDir pluginDir(QApplication::applicationDirPath());
 	QStringList plugin_filters; // filter for plugins, i.e "*.dll", "*.dylib"
 #if defined(Q_OS_WIN)
 	if (pluginDir.dirName().toLower() =="debug"
 		||pluginDir.dirName().toLower() =="release")
 		pluginDir.cdUp();
 	plugin_filters << "*.dll";
 	#elif defined(Q_OS_MAC)
 	if (pluginDir.dirName() =="MacOS"){
 		pluginDir.cdUp();
 		pluginDir.cdUp();
 		pluginDir.cdUp();
+	}
 	plugin_filters << "*.dylib";
 #endif
 	pluginDir.setNameFilters(plugin_filters);
 	if (!pluginDir.cd("models")) {
 		MyWarning::error("Directory 'models' not found!");
+	}
 	//QVector<SimPluginInterface *> plugins;
 	foreach (QString fileName, pluginDir.entryList(QDir::Files)){
 		QPluginLoader loader(pluginDir.absoluteFilePath(fileName));
 		if (SimPluginInterface *plugin =
 			qobject_cast<SimPluginInterface *>(loader.instance())) {
 			models.append(plugin);
 		} else {
 			MyWarning::warning("Could not load plugin %s",fileName.toStdString().c_str());
 			cerr << loader.errorString().toStdString().c_str() << endl;
 			MyWarning::warning("Could not load model %s: %s",fileName.toStdString().c_str(), loader.errorString().toStdString().c_str());
+		}
+	}
 	if (models.size()==0) {
 		MyWarning::error("No models could be loaded.");
+	}
+}
 void ModelCatalogue::LoadPlugin(const char *model) {
   QDir pluginDir(QApplication::applicationDirPath());
   QStringList plugin_filters; // filter for plugins, i.e "*.dll", "*.dylib"
 #if defined(Q_OS_WIN)
   if (pluginDir.dirName().toLower() =="debug"
       ||pluginDir.dirName().toLower() =="release")
     pluginDir.cdUp();
   //plugin_filters << "*.dll";
 #elif defined(Q_OS_MAC)
   if (pluginDir.dirName() =="MacOS"){
     pluginDir.cdUp();
     pluginDir.cdUp();
     pluginDir.cdUp();
+  }
   //plugin_filters << "*.dylib";
 #endif
   plugin_filters << model;
   pluginDir.setNameFilters(plugin_filters);
   if (!pluginDir.cd("models")) {
     MyWarning::error("Directory 'models' not found!");
+  }
   //QVector<SimPluginInterface *> plugins;
   QStringList modelnames=pluginDir.entryList(QDir::Files);
   if (modelnames.empty()) {
     MyWarning::error("Model %s not found - hint: do not include path in filename.",model);
+  }
   foreach (QString fileName, modelnames){
     QPluginLoader loader(pluginDir.absoluteFilePath(fileName));
     if (SimPluginInterface *plugin =
 	qobject_cast<SimPluginInterface *>(loader.instance())) {
       models.append(plugin);
       //MyWarning::warning("Successfully loaded model %s",fileName.toStdString().c_str());
     } else {
       MyWarning::warning("Could not load plugin %s",fileName.toStdString().c_str());
+    }
+  }
+}
 void ModelCatalogue::InstallFirstModel() {
 	InstallModel(models[0]);
+}
 void ModelCatalogue::PopulateModelMenu() {
 	foreach (SimPluginInterface *model, models) {
 		QAction *modelaction = new QAction(model->ModelID(), mainwin);
 		QVariant data;
 		data.setValue(model);
 		modelaction->setData(data);
 		mainwin->modelmenu->addAction(modelaction);
+	}
 	connect(mainwin->modelmenu, SIGNAL(triggered(QAction *)), this, SLOT(InstallModel(QAction *)) );
+}
 void ModelCatalogue::InstallModel(QAction *modelaction) {
 	QVariant data = modelaction->data();
 	SimPluginInterface *model = data.value<SimPluginInterface *>();
 	cerr << "You chose model " << model->ModelID().toStdString() << "!\n";
 	mesh->Clean();
 	InstallModel(model);
+}
 void ModelCatalogue::InstallModel(SimPluginInterface *plugin) {
   // make sure both main and plugin use the same static datamembers (ncells, nchems...)
   plugin->SetCellsStaticDatamembers(CellBase::GetStaticDataMemberPointer());
   mesh->SetSimPlugin(plugin);
   Cell::SetNChem(plugin->NChem());
   plugin->SetParameters(&par);
   if (mainwin) {
     mainwin->RefreshInfoBar();
     mainwin->Init(0);
+  }
   //	mesh->StandardInit();
+}

src/simplugin.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.
+ *
  */
 #ifndef _SIMPLUGIN_H_
 #define _SIMPLUGIN_H_
 #include <QtPlugin>
 #include <QMetaType>
 #include "cellbase.h"
 #include "wallbase.h"
 class Parameter;
 #include <QColor>
 #include <QString>
 class SimPluginInterface {
 public:
 	virtual QString ModelID(void) = 0;
 	virtual ~SimPluginInterface() { }
 	// Executed after the cellular mechanics steps have equillibrized
-	virtual void CellHouseKeeping(CellBase &c) = 0;
+	virtual void CellHouseKeeping(CellBase *c) = 0;
 	// Differential equations describing transport of chemicals from cell to cell
 	virtual void CelltoCellTransport(Wall *, double *dchem_c1, double *dchem_c2) = 0;
 	// Differential equations describing chemical reactions taking place at or near the cell walls
 	// (e.g. PIN accumulation)
 	virtual void WallDynamics(Wall *w, double *dw1, double *dw)  = 0;
 	// Differential equations describing chemical reactions inside the cells
 	virtual void CellDynamics(CellBase *c, double *dchem) = 0;
 	// to be executed after a cell division
-	virtual void OnDivide(ParentInfo &parent_info, CellBase &daughter1, CellBase &daughter2) = 0;
+	virtual void OnDivide(ParentInfo *parent_info, CellBase *daughter1, CellBase *daughter2) = 0;
 	// to be executed for coloring a cell
-	virtual void SetCellColor(CellBase &c, QColor &color) = 0;
+	virtual void SetCellColor(CellBase *c, QColor *color) = 0;
 	// Number of chemicals
 	virtual int NChem(void) = 0;
 	// For internal use; not to be redefined by end users
 	virtual void SetParameters(Parameter *pass_pars);// { par = pass_pars; }
 	virtual void SetCellsStaticDatamembers (CellsStaticDatamembers *cells_static_data_members_of_main);
 protected:
 	class Parameter *par;
 };
 Q_DECLARE_INTERFACE(SimPluginInterface,
-                    "nl.cwi.VirtualLeaf.SimPluginInterface/1.1")
+                    "nl.cwi.VirtualLeaf.SimPluginInterface/1.2")
 Q_DECLARE_METATYPE(SimPluginInterface *)
 #endif
@@ \ No newline at end of file @@

0 comments (0 inline, 0 general)