Files @ 6bcb69712a0e
Branch filter:

Location: EI/VirtualLeaf/src/nodeset.h

Roeland Merks
In response to referee's comment:

"However, (although it is probably not important for model developers), I was
still unable to load the model ‘Meinhardt 1976 with growth’ example from the
‘Models’ drop down menu, I got the ‘Fatal error’ message ‘stepwise underflow in

rkqs, with h=0.000000 and htry = 0.100000’. The model did work when I selected
the Meinhardt model in both the Models and the File -> Read leaf menus together,
it would be preferable if the models loaded from the Models menu worked

automatically. I am not sure that the update referred to in the author’s
response permits the loading of ‘My First Model in Virtual Leaf’ from the
‘Models’ drop down menu; I only got a cell that wiggled."

I made the following changes:

- meinhardt_plugin now reads "meinhardt_init.xml" by default
- changed the name of 'My First Model in Virtual Leaf’ to '0: Empty model template (does nothing)' to make it clear that the wiggle cell is the intended behavior for this model example.
- Added default parameter files for Tutorial1A-D and Tutorial2 to prevent variable behavior depending on the last parameter settings used.


--
user: Roeland Merks <roeland.merks@cwi.nl>
branch 'default'
added data/leaves/tutorial1_init.xml
added data/leaves/tutorial2_init.xml
changed data/leaves/auxin_growth.xml
changed data/leaves/meinhardt_init.xml
changed src/TutorialCode/Tutorial0/tutorial0.cpp
changed src/TutorialCode/Tutorial1A/tutorial1A.h
changed src/TutorialCode/Tutorial1B/tutorial1B.h
changed src/TutorialCode/Tutorial1C/tutorial1C.h
changed src/TutorialCode/Tutorial1D/tutorial1D.h
changed src/TutorialCode/Tutorial2/tutorial2.h
changed src/build_models/meinhardtplugin.h
/*
 *
 *  $Id$
 *
 *  This file is part of the Virtual Leaf.
 *
 *  VirtualLeaf 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.
 *
 *  VirtualLeaf 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 _NODESET_H_
#define _NODESET_H_

#include <algorithm>
#include <numeric>
#include <list>
#include <iterator>
#include "node.h"

class NodeSet : public list<Node *> {

 public:
  NodeSet(void) {
    done = false;
  }

  inline bool DoneP(void) { return done; }
  inline void ResetDone(void) { done = false; }

  void AddNode(Node * n) {
    push_back(n);
    n->node_set = this;
  }

  list <Cell *> getCells(void) {

    list<Cell *> cellset;

    for (list<Node *>::const_iterator i=begin();
	 i!=end();
	 i++) {
      transform ( (*i)->owners.begin(), (*i)->owners.end(), back_inserter( cellset ) , mem_fun_ref ( &Neighbor::getCell ));
    }

    cellset.sort();

    // remove all non-unique elements
    // (unique moves all unique elements to the front and returns an iterator to the end of the unique elements;
    // so we simply erase all remaining elements.
    cellset.erase(::unique(cellset.begin(), cellset.end() ), cellset.end() );

    // remove boundary_polygon
    cellset.erase( find_if ( cellset.begin(), cellset.end(), mem_fun( &Cell::BoundaryPolP  ) ) );
    return cellset;
  }

  void CleanUp(void) {

    // remove double Nodes from the set
    sort();
    erase(::unique(begin(), end() ), 
	  end() );


  }
  void print(ostream &os) const {
    transform( begin(), end(), ostream_iterator<int>( os, " " ), mem_fun( &Node::Index ) );
  }

  /*! Attempt a move over (rx, ry)
    reject if energetically unfavourable.
  */

  void AttemptMove(double rx, double ry) {

    done = true;
    // 1. Collect list of all attached to the nodes in the set
    list<Cell *> celllist = getCells();

    // 2. Sum the current energy of these cells
    double old_energy=0.;
    double sum_stiff = 0.;
    for ( list<Cell *>::const_iterator i = celllist.begin(); 
	  i!=celllist.end();
	  ++i ) {

      old_energy += (*i)->Energy();
      sum_stiff += (*i)->Stiffness();
    }

    // 3. (Temporarily) move the set's nodes.
    for ( list<Node *>::iterator n = begin();
	  n!=end();
	  ++n ) {

      (*n)->x+=rx;
      (*n)->y+=ry;
      // (*n)->z += rz;
    }

    // 4. Recalculate the energy
    double new_energy = 0.;
    for ( list<Cell *>::const_iterator i = celllist.begin(); 
	  i!=celllist.end();
	  ++i ) {

      new_energy += (*i)->Energy();
    }


    // 5. Accept or Reject DeltaH
    double dh = new_energy - old_energy;

    list<int> new_areas;

    // cerr << "Nodeset says: dh = " << dh << " ...";
    if (dh < -sum_stiff || RANDOM()<exp((-dh - sum_stiff)/par.T) ) {

      // ACCEPT
      // recalculate areas of cells
      for_each ( celllist.begin(), celllist.end(), mem_fun ( &Cell::RecalcArea ) );

      // recalculate integrals of cells
      for_each ( celllist.begin(), celllist.end(), mem_fun ( &Cell::SetIntegrals ) );


    } else {
      // REJECT

      // Move the set's nodes back to their original position
      for ( list<Node *>::iterator n = begin();
	    n!= end();
	    ++n ) {

	(*n)->x-=rx;
	(*n)->y-=ry;
      }
    }
  }

  void XMLAdd(xmlNode *root) const;
  void XMLRead(xmlNode *root, Mesh *m);
 private:
  bool done;
};

ostream &operator<<(ostream &os, const NodeSet &ns);
#endif

/* finis */