diff --git a/src/cell.cpp b/src/cell.cpp new file mode 100644 --- /dev/null +++ b/src/cell.cpp @@ -0,0 +1,1989 @@ +/* + * + * 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 . + * + * Copyright 2010 Roeland Merks. + * + */ + +#include +#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 + +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::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::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::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 superfluous_nodes; + for (list::iterator n=nodes.begin(); + n!=nodes.end(); + n++) { + + Node &no(*(*n)); + // locate myself in the node's owner list + list::iterator cellpos; + bool cell_found=false; + for (list::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::iterator n=nodes.begin(); + n!=nodes.end(); + n++) { + + Node &no(*(*n)); + + if (no.DeadP()) { + + list::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::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::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::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::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::iterator previous_iterator=i; + previous_iterator--; + previous=*previous_iterator; + } else { + previous=nodes.back(); + } + + Node *next; + list::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::iterator i=nodes.begin(); + i!=nodes.end(); + i++) { + + Vector v3 = *(*i); + list::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( 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;ichem[i]=chem[i]; + } + + //extern double auxin_account; + //auxin_account += daughter->chem[0]; + + for (int i=0;inew_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() <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::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::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 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 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 edge_owners; + list::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::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::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::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 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::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::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::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::iterator i=start; + while ( i!=stop) { + + // give the node to the daughter + // (find references to parent cell from this node, + // and remove them) + list::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 new_nodes_parent; + + // half of the nodes stay with the parent + { + list::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::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::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 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 copy_walls = walls; + for (list::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::iterator i=broken_neighbors.begin(); + i!=broken_neighbors.end();i++) { + ((Cell *)(*i))->ConstructNeighborList(); + } + + + ConstructNeighborList(); + daughter->ConstructNeighborList(); + + 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::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 > length_edges; + vector 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::const_iterator i=nodes.begin(); + i!=nodes.end(); + i++) { + + //if ((*i)->Fixed()) return; // commented out 01/12/05 + for (list::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 (*i, n->nb1) ); + length_edges.push_back( pair (*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::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::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()::const_iterator nb_it = neighbors.begin(); + for (vector::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 area1; + vector area2; + m->ExtractFromCells( mem_fun_ref(&Cell::Area), back_inserter(area1) ); + m->ExtractFromCells( mem_fun_ref(&Cell::CalcArea), back_inserter(area2)); + vector::iterator i=area1.begin(); + vector::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::const_iterator i=nodes.begin(); + i!=nodes.end(); + i++) { + + for (list::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::const_iterator i=nodes.begin(); + i!=nodes.end(); + i++) { + + list::const_iterator j=i; + ++j; + for (; + j!=nodes.end(); + j++) + { + + Vector v1 = *(*i); + list::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::const_iterator moving_node_ind_pos = find (nodes.begin(),nodes.end(),moving_node_ind);\n"; + list::const_iterator moving_node_ind_pos = find (nodes.begin(),nodes.end(),moving_node_ind); + + list::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::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::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::const_iterator i=nodes.begin(); + i!=nodes.end(); + i++) + { + + Vector v3 = *(*i); + list::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 corner_points; + + for (list::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::const_iterator nb = (--corner_points.end()); + + // loop over list, + for (list::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 owning_cells; + Node &n(*(*i)); + + for (list::const_iterator j=n.owners.begin(); + j!=n.owners.end(); + j++) { + owning_cells.push_back(j->cell); + } + + Node &n2(*(*nb)); + for (list::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 duplicates; + list::const_iterator prevj = (--owning_cells.end()); + for (list::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::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::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::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;clength * ( 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::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); + + 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::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::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( 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::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::const_iterator first_node_edge = find(nodes.begin(), nodes.end(), (*de)->n1); + list::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::const_iterator n=++first_node_edge; + n!=second_node_edge_plus_1; + ++n ) { + + if (n==nodes.end()) n=nodes.begin(); /* wrap around */ + + + list::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::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); + +}