diff --git a/src/cell.cpp b/src/cell.cpp --- a/src/cell.cpp +++ b/src/cell.cpp @@ -40,1510 +40,1461 @@ extern Parameter par; double Cell::factor=1.; double Cell::offset[3]={0,0,0}; -Cell::Cell(void) : CellBase() { - - m=0; - +Cell::Cell(void) : CellBase() +{ + m=0; } -Cell::Cell(double x, double y, double z) : CellBase(x,y,z) { - - 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; +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 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); - - } - prev_cross_z=cross.z; - } - - DivideWalls(new_node_locations, centroid, centroid+axis); - +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); + + } + prev_cross_z=cross.z; + } + + DivideWalls(new_node_locations, centroid, centroid+axis); } -double Cell::MeanArea(void) { - return m->MeanArea(); + +double Cell::MeanArea(void) +{ + return m->MeanArea(); } -void Cell::Apoptose(void) { - - // First kill walls - #ifdef QDEBUG - qDebug() << "This is cell " << Index() << endl; - qDebug() << "Number of walls: " << walls.size() << endl; - #endif - for (list::iterator w=walls.begin(); w!=walls.end(); w++) { - #ifdef QDEBUG - qDebug() << "Before apoptosis, wall " << (*w)->Index() << " says: c1 = " - << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl; - #endif - } - 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; - } - - #ifdef QDEBUG - if (illegal_flag && (*w)->c1==(*w)->c2) { - qDebug() << "I created an illegal wall." << endl; - } - #endif +void Cell::Apoptose(void) +{ + // First kill walls +#ifdef QDEBUG + qDebug() << "This is cell " << Index() << endl; + qDebug() << "Number of walls: " << walls.size() << endl; +#endif + for (list::iterator w=walls.begin(); w!=walls.end(); w++) { +#ifdef QDEBUG + qDebug() << "Before apoptosis, wall " << (*w)->Index() << " says: c1 = " + << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl; +#endif + } + 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) { - if ( ((*w)->N1()->DeadP() || (*w)->N2()->DeadP()) || - ((*w)->C1() == (*w)->C2() ) ){ - // kill wall - #ifdef QDEBUG - qDebug() << "Killing wall." << endl; - #endif - (*w)->Kill(); + // 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; + } + +#ifdef QDEBUG + if (illegal_flag && (*w)->c1==(*w)->c2) { + qDebug() << "I created an illegal wall." << endl; + } +#endif + + if ( ((*w)->N1()->DeadP() || (*w)->N2()->DeadP()) || + ((*w)->C1() == (*w)->C2() ) ){ + // kill wall +#ifdef QDEBUG + qDebug() << "Killing wall." << endl; +#endif + (*w)->Kill(); - #ifdef QDEBUG - if ((*w)) { - qDebug() << "Wall " << (*w)->Index() << " says: c1 = " - << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl; - } - #endif - (*w)=0; - } else { - #ifdef QDEBUG - qDebug() << "Not killing wall." << endl; - qDebug() << "Wall " << (*w)->Index() << " says: c1 = " - << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl; - #endif - } - - } - 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); - } - } - } - - - - // mark cell as dead - MarkDead(); +#ifdef QDEBUG + if ((*w)) { + qDebug() << "Wall " << (*w)->Index() << " says: c1 = " + << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl; + } +#endif + (*w)=0; + } else { +#ifdef QDEBUG + qDebug() << "Not killing wall." << endl; + qDebug() << "Wall " << (*w)->Index() << " says: c1 = " + << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl; +#endif + } + } + 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); + } + } + } + + // 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; - // } +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; + } + (*i)->owners.push_back( Neighbor( this, previous, next ) ); + } } +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; + +#ifdef QDEBUG + qDebug() << "Cell " << Index() << " is doing DivideOverGivenLine" << endl; +#endif + 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; + -/*! \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; - - #ifdef QDEBUG - qDebug() << "Cell " << Index() << " is doing DivideOverGivenLine" << endl; - #endif - 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; - - - //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); - - } - } - - #ifdef QDEBUG - if (new_node_locations.size()<2) { - qDebug() << "Line does not intersect with two edges of Cell " << Index() << endl; - qDebug() << "new_node_locations.size() = " << new_node_locations.size() << endl; - return false; - } - - ItList::iterator i = new_node_locations.begin(); - list< Node *>::iterator j; - qDebug() << "-------------------------------" << endl; - qDebug() << "Location of new nodes: " << (**i)->Index() << " and "; + //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); + + } +} + +#ifdef QDEBUG + if (new_node_locations.size()<2) { + qDebug() << "Line does not intersect with two edges of Cell " << Index() << endl; + qDebug() << "new_node_locations.size() = " << new_node_locations.size() << endl; + return false; + } - ++i; - j = *i; - if (j==nodes.begin()) j=nodes.end(); j--; - - qDebug() << (*j)->Index() << endl; - qDebug() << "-------------------------------" << endl; - - if ( **new_node_locations.begin() == *j ) { - qDebug() << "Rejecting proposed division (cutting off zero area)." << endl; - return false; - } - #endif +ItList::iterator i = new_node_locations.begin(); +list< Node *>::iterator j; +qDebug() << "-------------------------------" << endl; +qDebug() << "Location of new nodes: " << (**i)->Index() << " and "; + +++i; +j = *i; +if (j==nodes.begin()) j=nodes.end(); j--; - DivideWalls(new_node_locations, v1, v2, fix_cellwall, node_set); - - return true; - +qDebug() << (*j)->Index() << endl; +qDebug() << "-------------------------------" << endl; + +if ( **new_node_locations.begin() == *j ) { + qDebug() << "Rejecting proposed division (cutting off zero area)." << endl; + return false; + } +#endif + +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]; - } - - daughter->cell_type = cell_type; - - 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; - - - // 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) - - #ifdef QDEBUG - qDebug() << "new_node_locations.size() = " << new_node_locations.size() <index = " << daughter->index << endl; - qDebug() << "cells.size() = " << m->cells.size() << endl; - #endif +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]; + } + + daughter->cell_type = cell_type; + + 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; + + // 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) + +#ifdef QDEBUG + qDebug() << "new_node_locations.size() = " << new_node_locations.size() <index = " << daughter->index << endl; + qDebug() << "cells.size() = " << m->cells.size() << endl; +#endif - 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; + 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; + 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 << endl ; - } 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; + + 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 << endl ; + } 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] ); } - - - 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 - } - - // if new node is inserted into the boundary - // it will be part of the boundary, too + + // 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 + } + + // 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; - 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. + // 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) { - #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(); + m->boundary_polygon->nodes.insert(((ins_pos--)!=m->boundary_polygon->nodes.begin()?ins_pos:(--m->boundary_polygon->nodes.end())), new_node_ind[i]); - // 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 + // .. 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; - 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 ) ); + // 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 +#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); + } + } - // 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 - 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); - } - } + // 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()) { - #ifdef QDEBUG - qDebug() << "The edge owners list has: " << edge_owners.size() << " elements" << endl; - #endif + // 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(); - // 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(); + // 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; - // (can we have more than one neighboring cell here??) - if (c!=owners.end()) { - neighbor_cell = c->cell; - if (!c->cell->BoundaryPolP()) { + // 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 ) ); - // 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()) { +#ifdef QDEBUG + qDebug() << "Whoops, wall element not found...!" << endl; + qDebug() << "Cell ID: " << neighbor_cell->Index() << endl; + qDebug() << "My cell ID: " << Index() << endl; +#endif + } 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; - 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()) { - #ifdef QDEBUG - qDebug() << "Whoops, wall element not found...!" << endl; - qDebug() << "Cell ID: " << neighbor_cell->Index() << endl; - qDebug() << "My cell ID: " << Index() << endl; - #endif - } 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]; - } - - - // 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()) { - #ifdef QDEBUG - qDebug() << "not found" << endl; - #endif - abort(); - } - - (*i)->owners.erase(neighb_with_this_cell); - - daughter->nodes.push_back( *i ); - - - i++; - if (i==nodes.end()) - i=nodes.begin(); - }; + // 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]; + } + + + // 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; + } } - - // 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(); - }; - } + } + } + } // 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()) { + +#ifdef QDEBUG + qDebug() << "not found" << endl; +#endif + abort(); + } + + (*i)->owners.erase(neighb_with_this_cell); + + daughter->nodes.push_back( *i ); + - // 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] ); - } - + 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 { - node->fixed = 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 ); - } - - } + // 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 { + node->fixed = 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]); - } - } - - //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); - - daughter->div_counter=(++div_counter); - - + 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]); + } + } + + //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); + + 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; - } + + 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++) { - +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++) { + + + for (list::const_iterator n=(*i)->owners.begin(); + n!=(*i)->owners.end(); + n++) { + + if (n->getCell()!=this) { + 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()); + } + } + } - for (list::const_iterator n=(*i)->owners.begin(); - n!=(*i)->owners.end(); - n++) { - - if (n->getCell()!=this) { - 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; - - } else { - - Move ( -1*movement); - - } - - return dh; + // 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; + + } 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); +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; +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/ +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++) + { - 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; - - - if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) { - //cerr << "ua = " << ua << ", ub = " << ub << endl; - return true; - } - } - } - - return false; + 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; + + + 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; - - 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; + + 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++) +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; - - if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) { - return true; - } + 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; + + if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) { + return true; + } } - - return false; - - + 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) { - + + 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; - 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 )); + } + } + + + // One Wall for each neighbor, so we should be able to correctly construct neighbor lists here. + if (duplicates[0]==this) { + AddWall( new Wall(*nb,*i,duplicates[0],duplicates[1]) ); + if (!duplicates[1]->BoundaryPolP()) { - // 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) { - AddWall( new Wall(*nb,*i,duplicates[0],duplicates[1]) ); - if (!duplicates[1]->BoundaryPolP()) { - - neighbors.push_back(duplicates[1]); - } - } else { - AddWall ( new Wall(*nb,*i,duplicates[1],duplicates[0]) ); - if (!duplicates[0]->BoundaryPolP()) { - neighbors.push_back(duplicates[0]); - - } - } - } - + neighbors.push_back(duplicates[1]); + } + } else { + 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 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] ); - - #ifdef QDEBUG - if ((*i)->c1!=this) { - qDebug() << "Warning, bad cells boundary: " << (*i)->c1->Index() << ", " << index << endl; - } - #endif - - flux[c] += phi; - } - } - +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] ); + +#ifdef QDEBUG + if ((*i)->c1!=this) { + qDebug() << "Warning, bad cells boundary: " << (*i)->c1->Index() << ", " << index << endl; + } +#endif + + flux[c] += phi; + } + } } @@ -1552,46 +1503,46 @@ void Cell::Flux(double *flux, double *D) #include "canvas.h" -void Cell::Draw(QGraphicsScene *c, QString tooltip) { - - // Draw the cell on a QCanvas object - - if (DeadP()) { - #ifdef QDEBUG - qDebug() << "Cell " << index << " not drawn, because dead." << endl; - #endif - 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::Draw(QGraphicsScene *c, QString tooltip) +{ + + // Draw the cell on a QCanvas object + + if (DeadP()) { +#ifdef QDEBUG + qDebug() << "Cell " << index << " not drawn, because dead." << endl; +#endif + 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(); } @@ -1599,234 +1550,227 @@ void Cell::DrawCenter(QGraphicsScene *c) // Maginfication derived similarly to that in nodeitem.cpp // Why not use Cell::Magnification()? const double mag = par.node_mag; - - // construct an ellipse + + // 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); + 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; - + + for (list::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(((offset[0]+i->x)*factor), - ((offset[1]+i->y)*factor) ); - } - + NodeItem *item = new NodeItem ( &(*i), c ); + item->setColor(); + item->setZValue(5); + item->show(); + item ->setPos(((offset[0]+i->x)*factor), + ((offset[1]+i->y)*factor) ); + } } void Cell::DrawIndex(QGraphicsScene *c) const { - - DrawText( c, QString("%1").arg(index)); + + 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->show(); - ctext ->setPos(((offset[0]+centroid.x)*factor), - ((offset[1]+centroid.y)*factor) ); - + + 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->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 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; + 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(); - + 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 { c = NULL; // assignment merely to obviate compilation warning - MyWarning::warning("Sorry, Cell::DrawStrain temporarily not implemented."); + MyWarning::warning("Sorry, Cell::DrawStrain temporarily not implemented."); } -void Cell::DrawFluxes(QGraphicsScene *c, double arrowsize) { - - // get the mean flux through this cell - 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::DrawFluxes(QGraphicsScene *c, double arrowsize) +{ + + // get the mean flux through this cell + 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 ) ); + + 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()) ); - + + DrawText(c, QString("%1").arg(walls.size()) ); } -#endif +#endif // QTGRAPHICS ! /*! \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 - } + + 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 - #ifdef QDEBUG - if (w->c1 == w->c2 ){ - qDebug() << "Wall between identical cells: " << w->c1->Index()<< endl; - } - #endif +void Cell::AddWall( Wall *w ) +{ + + // if necessary, we could try later inserting it at the correct position +#ifdef QDEBUG + if (w->c1 == w->c2 ){ + qDebug() << "Wall between identical cells: " << w->c1->Index()<< endl; + } +#endif - // 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); - } - + // 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 ) - ); - +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) +{ -void Cell::EmitValues(double t) { - - // cerr << "Attempting to emit " << t << ", " << chem[0] << ", " << chem[1] << endl; - emit ChemMonValue(t, chem); - + // cerr << "Attempting to emit " << t << ", " << chem[0] << ", " << chem[1] << endl; + emit ChemMonValue(t, chem); } + +/* finis */