CellDivider.cpp

Go to the documentation of this file.

/*
 * Copyright 2011-2016 Universiteit Antwerpen
 *
 * Licensed under the EUPL, Version 1.1 or  as soon they will be approved by
 * the European Commission - subsequent versions of the EUPL (the "Licence");
 * You may not use this work except in compliance with the Licence.
 * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the Licence is distributed on an "AS IS" basis,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the Licence for the specific language governing
 * permissions and limitations under the Licence.
 */
#include "CellDivider.h"

#include "bio/Cell.h"
#include "bio/CellAttributes.h"
#include "bio/Edge.h"
#include "bio/Mesh.h"
#include "bio/NeighborNodes.h"
#include "bio/Node.h"
#include "bio/Wall.h"
#include "math/array3.h"
#include "math/constants.h"
#include "math/Geom.h"
#include "math/RandomEngine.h"
#include "model/ComponentFactoryProxy.h"
#include "sim/CoreData.h"
#include "util/container/CircularIterator.h"
#include "util/container/ConstCircularIterator.h"
#include "util/misc/Exception.h"
#include "util/misc/log_debug.h"

#include <trng/uniform01_dist.hpp>
#include <array>
#include <cassert>
#include <list>
#include <memory>
#include <string>
#include <vector>

using namespace std;
using namespace SimPT_Sim;
using namespace boost::property_tree;
using namespace SimPT_Sim::Container;

namespace {

void AssignExistingNodesToSplitCells(Mesh* mesh, Cell* this_cell,
        vector<Node*>& parent_nodes, vector<Node*>& daughter_nodes,
        const array<Node*, 2>& new_node_index)
{
        const auto start = find(parent_nodes.begin(), parent_nodes.end(), new_node_index[0]);
        assert(start != parent_nodes.end()
                && "The new intersection should already be in the list of parent nodes.");

        rotate(parent_nodes.begin(), start, parent_nodes.end());

        const auto stop = find(parent_nodes.begin(), parent_nodes.end(), new_node_index[1]);
        assert(stop != parent_nodes.end()
                && "The new intersection should already be in the list of parent nodes.");
        assert(stop != parent_nodes.begin()
                && "The two intersections are the same.");

        auto is_nb = [this_cell](const NeighborNodes& n) {
                return (n.GetCell()->GetIndex() == this_cell->GetIndex());
        };

        for (auto it = parent_nodes.begin() + 1; it != stop; it++) {
                auto& owners = mesh->GetNodeOwningNeighbors(*it);
                auto nb_it = find_if(owners.begin(), owners.end(), is_nb);
                assert((nb_it != owners.end()) && "Node manipulation failure!");
                owners.erase(nb_it);
                daughter_nodes.push_back(*it);
        }
        parent_nodes.erase(parent_nodes.begin(), stop + 1);
}

vector<Node*> ComputeIntersectNodes(Cell* cell, const array<double, 3>& axis)
        {
        const array<double, 3> centroid   = cell->GetCentroid();
        const array<double, 3> cross_prod = CrossProduct(axis, centroid - *(cell->GetNodes().back()));
        double prev_cross_z = cross_prod[2];
        vector<Node*> intersect_nodes;
        vector<Node*>& nodes = cell->GetNodes();
        for (unsigned int i = 0U; i < nodes.size(); i++) {
                // cross product to detect position of division
                const array<double, 3> cross = CrossProduct(axis, centroid - *nodes.at(i));
                // If the cross product changes sign, the line from centroid to node
                // has crossed the axis going from the previous node to the current one.
                if (cross[2] * prev_cross_z < 0) {
                        intersect_nodes.push_back(nodes.at(i));
                }
                prev_cross_z = cross[2];
        }
        return intersect_nodes;
        }

vector<Node*> ComputeTwoIntersectNodes(Cell* cell, const array<double, 3>& axis, std::vector<Node*> intersect_nodes)
        {
        std::vector<Node*> new_intersect_nodes;

        // ------------------------------------------------------------------------
        // calculate intersection points of axis
        // & calculate on which side of the long axis each intersection point lies
        // ------------------------------------------------------------------------
        std::vector<Node*>& nodes = cell->GetNodes();
        GeoData geo_dat = cell->GetGeoData();
        std::array<double, 3> long_axis = std::get<1>(geo_dat.GetEllipseAxes());

        const auto from = cell->GetCentroid();
        const auto to   = from + axis;

        std::map<Node*, double> intersect_points;
        for (Node* nn_ptr : intersect_nodes) {
                // -----------------------------------------------------------------------
                // Intersection between division axis and edge from this node
                // to its predecessor.
                // -----------------------------------------------------------------------
                const auto nn_it    = make_const_circular(nodes, find(nodes.begin(), nodes.end(), nn_ptr));
                const auto pr_it    = prev(nn_it);
                const auto n_vec    = Geom::Intersection(from, to, *(*nn_it), *(*pr_it));

                // ------------------------------
                // Calculate distance to centroid
                // ------------------------------
                double dist = sqrt(pow((cell->GetCentroid()[0] - n_vec[0]), 2) +
                        pow((cell->GetCentroid()[1] - n_vec[1]), 2));

                // ------------------------------------------------------------
                // Calculate which side of the long axis the intersection is on
                // ------------------------------------------------------------
                const array<double, 3> cross = CrossProduct(long_axis, cell->GetCentroid() - *nn_ptr);
                int sign = SimPT_Sim::Util::signum(cross[2]);

                // ----------------------------------------------------------------
                // Set distance - if sign or cross was -, + if sign of cross was +
                // ----------------------------------------------------------------
                dist *= sign;
                intersect_points[nn_ptr] = dist;
        }

        // --------------------------------------------------------------------
        // Choose closest intersection points on either side of the long axis
        // --------------------------------------------------------------------
        Node* node1;
        double node1_dist = std::numeric_limits<double>::infinity();
        Node* node2;
        double node2_dist = std::numeric_limits<double>::infinity();
        for (auto entry : intersect_points) {
                if(entry.second >= 0) {
                        if(entry.second < node1_dist) {
                                node1 = entry.first;
                                node1_dist = entry.second;
                        }
                } else {
                        if(abs(entry.second) < node2_dist) {
                                node2 = entry.first;
                                node2_dist = abs(entry.second);
                        }
                }
        }

        new_intersect_nodes.push_back(node1);
        new_intersect_nodes.push_back(node2);
        return new_intersect_nodes;

        }

} // anonymous namespace


namespace SimPT_Sim {

CellDivider::CellDivider(const CoreData& cd)
{
        Initialize(cd);
}

CellDivider::CellDivider(Mesh* mesh)
{
        Initialize(mesh);
}

void CellDivider::AddSharedWallNodes(array<array<double, 3>, 2>& new_node_vector,
        array<Node*, 2>& new_node_index, vector<Node*>& parent_nodes, vector<Node*>& daughter_nodes)
{
        //--------------------------------------------------------------------------
        // Insert the nodes that start the shared wall, but remember these are different
        // for parent and daughter because their traversal order is different in parent and daughter.
        //--------------------------------------------------------------------------
        parent_nodes.push_back(new_node_index[0]);
        daughter_nodes.push_back(new_node_index[1]);

        //--------------------------------------------------------------------------
        // Skip this for some of the models.
        //--------------------------------------------------------------------------
        if ( m_model_name != "Wortel"
                && m_model_name != "Maize"
                        && m_model_name != "MaizeGRN"
                                && m_model_name != "WortelLight"
                                        && m_model_name != "Blad"
                                                && m_model_name != "WrapperModel") {

                // Add intermediate, extra nodes along the shared wall. Calculate the
                // distance between new nodes to estimate the number of intermediate nodes.
                // The factor 4.0 is there to keep tension on the walls;
                // this is a hidden parameter and should be made explicit later on.

                const double dist      = Norm(new_node_vector[1] - new_node_vector[0]);
                const double tnd       = abs(m_target_node_distance);
                const int n            = (tnd > SimPT_Sim::Util::small_real()) ? static_cast<int>((dist/tnd) / 4.0 + 0.5) : 0;
                const double new_dist  = dist / static_cast<double>(n + 1);
                const auto nodevec     = Normalize(new_node_vector[1] - new_node_vector[0]);

                // Note that wall nodes need to run in inverse order in parent
                vector<Node*>::iterator ins_pos = daughter_nodes.end();
                for (int i = 1; i <= n; i++) {
                        const array<double, 3> vec = new_node_vector[0] + i * new_dist * nodevec;
                        auto node = m_mesh->BuildNode(vec, false);
                        ins_pos = daughter_nodes.insert(ins_pos, node);
                        parent_nodes.push_back(node);
                }
        }

        //--------------------------------------------------------------------------
        // Insert the final nodes of the shared wall, again taking traversal order into account.
        //--------------------------------------------------------------------------
        daughter_nodes.push_back(new_node_index[0]);
        parent_nodes.push_back(new_node_index[1]);
}

void CellDivider::AddSplitWalls(Cell* this_cell, Cell* neighbor_cell, int i,
        array<Node*, 2>& new_node_ptr, array<Wall*, 4>& div_wall,
        array<double, 2>& orig_length, array<double, 2>& orig_rest_length, array<double, 2>& orig_rest_length_init)
{
        auto this_cell_walls = this_cell->GetWalls();
        assert((neighbor_cell != nullptr) && "There has to be a neighbor_cell!");

        if (neighbor_cell && !this_cell_walls.empty() /* && !neighbor_cell->BoundaryPolP() */) {

                const auto is_of_nb =[neighbor_cell](Wall* w) {
                        return (w->GetC1() == neighbor_cell) || (w->GetC2() == neighbor_cell);
                };

                //--------------------------------------------------------------------------
                // Find the correct wall element. This is the wall between this cell and neighbor
                // cell containing the intersection node.
                // All intersection nodes have already been added to this cell.
                //--------------------------------------------------------------------------
                auto w = this_cell_walls.begin();
                bool ok = false;
                do {
                        w = find_if(w, this_cell_walls.end(), is_of_nb);

                        if (w == this_cell_walls.end()) {
                                break;
                        }

                        auto nodes = (*w)->GetC1()->GetNodes();
                        rotate(nodes.begin(), find(nodes.begin(), nodes.end(), (*w)->GetN1()), nodes.end());
                        auto secondNode = find(nodes.begin(), nodes.end(), (*w)->GetN2());

                        for (auto it = nodes.begin(); it != next(secondNode); it++) {
                                if (*it == new_node_ptr[i]) {
                                        ok = true;
                                        break;
                                }
                        }
                } while (!ok && ++w != this_cell_walls.end());  //Short circuiting to have the right w if ok.
                assert((w != this_cell_walls.end()) && "No split wall element found!");

                //--------------------------------------------------------------------------
                // 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_ptr[i] != (*w)->GetN1()) && (new_node_ptr[i] != (*w)->GetN2())) {
                        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)->SetLengthDirty();
                        orig_length[i] = (*w)->GetLength();
                        orig_rest_length[i] = (*w)->GetRestLength();
                        orig_rest_length_init[i] = (*w)->GetRestLengthInit();
                        if ((*w)->GetC1() == this_cell) {
                                new_wall = m_mesh->BuildWall((*w)->GetN1(),
                                        new_node_ptr[i], this_cell, neighbor_cell);
                        } else {
                                new_wall = m_mesh->BuildWall((*w)->GetN1(),
                                        new_node_ptr[i], neighbor_cell, this_cell);
                        }
                        (*w)->SetN1(new_node_ptr[i]);
                        (*w)->SetLengthDirty();
                        (*w)->SetRestLength((orig_rest_length[i] / orig_length[i])*(*w)->GetLength());
                        (*w)->SetRestLengthInit((orig_rest_length_init[i] / orig_length[i])*(*w)->GetLength());

                        //--------------------------------------------------------------------------
                        //test SETWALLSTRENGTH transplant
                        //--------------------------------------------------------------------------
                        (*w)->SetStrength((*w)->GetStrength());
                        new_wall->SetStrength((*w)->GetStrength());
                        assert(new_wall->GetN1() != new_wall->GetN2() && "What the hell!");

                        //--------------------------------------------------------------------------
                        // Give wall elements to appropriate cells.
                        //--------------------------------------------------------------------------
                        if (m_copy_wall) {
                                new_wall->WallAttributes::operator=(**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
                                //--------------------------------------------------------------------------
                                assert((m_cd.m_random_engine != nullptr) && "Random engine not initialized.");
                                if (m_uniform_generator() < 0.5) {
                                        new_wall->WallAttributes::Swap(**w);
                                }
                        }
                        //--------------------------------------------------------------------------
                        // Remember the addresses of the new walls.
                        //--------------------------------------------------------------------------
                        this_cell->AddWall(new_wall);
                        neighbor_cell->AddWall(new_wall);
                        div_wall[2 * i + 0] = *w;
                        div_wall[2 * i + 1] = new_wall;

                        //--------------------------------------------------------------------------
                        // Correct the walls associated with the nodes on the new wall
                        // Old wall should also be updated, as a node might have been added to split the wall
                        //--------------------------------------------------------------------------
                        m_mesh->UpdateNodeOwningWalls(*w);
                        m_mesh->UpdateNodeOwningWalls(new_wall);
                }
        }
}

void CellDivider::ComputeIntersection(Cell* cell, vector<Node*> intersect_nodes,
        const array<double, 3>& from, const array<double, 3>& to,
        array<array<double, 3>, 2>& new_node_vector, array<Node*, 2>& new_node_ptr,
        array<int, 2>& new_node_flag, array<Edge, 2>& intersect_edge)
{
        int      nnc = 0;

        //--------------------------------------------------------------------------
        // Determine the intersection points and find out whether they (almost)
        // coincide with existing points in which we will use those existing
        // nodes to build the shared wall.
        //--------------------------------------------------------------------------
        vector<Node*>& nodes = cell->GetNodes();
        for (Node* nn_ptr : intersect_nodes) {

                //--------------------------------------------------------------------------
                // Intersection between division axis and edge from this node
                // to its predecessor. Remember the edge gets divided.
                //--------------------------------------------------------------------------
                const auto nn_it    = make_const_circular(nodes, find(nodes.begin(), nodes.end(), nn_ptr));
                const auto pr_it    = prev(nn_it);
                const auto n_vec    = Geom::Intersection(from, to, *(*nn_it), *(*pr_it));
                intersect_edge[nnc] = Edge(*pr_it, *nn_it);

                //--------------------------------------------------------------------------
                // Insert this new node if it is far enough (5% of element length)
                // from one of the two existing nodes, else use existing node.
                // new_node_flag == 0 : an actual new node needs to be created
                // new_node_flag == 1 : node pointed to by new_node_location iterator is used
                // new_node_flag == 2 : predecessor of node pointed to by new_node_location iterator is used
                //--------------------------------------------------------------------------
                const array<double, 3> nn = *(*nn_it);
                const array<double, 3> pr = *(*pr_it);
                const double tol = m_collapse_node_threshold * Norm(nn - pr);
                if (Norm(nn - n_vec) <= tol) {
                        new_node_flag[nnc]    = 1;
                        new_node_vector[nnc]  = nn;
                        new_node_ptr[nnc]     = *nn_it;
                } else if (Norm(pr - n_vec) <= tol) {
                        new_node_flag[nnc]    = 2;
                        new_node_vector[nnc]  = pr;
                        new_node_ptr[nnc]     = *pr_it;
                } else {
                        new_node_flag[nnc]    = 0;
                        new_node_vector[nnc]  = n_vec;
                        new_node_ptr[nnc]     = nullptr;
                }

                nnc++;
        }
}

unsigned int CellDivider::DivideCells()
{
        unsigned int division_count = 0U;

        // Copy current cells (vector of current cells gets updated
        // by division so you cannot loop over it).
        const vector<Cell*> current_cells { m_mesh->GetCells() };
        // CellSplit tuple indicates whether to divide, whether
        // to use a fixed axis, and what that axis must be.
        for (auto const& cell : current_cells) {
                auto tup = m_cell_split(cell);
                if (get<0>(tup)) {
                        if (get<1>(tup)) {
                                DivideOverAxis(cell, get<2>(tup));
                        } else {
                                const auto axis = Orthogonalize(get<1>(cell->GetGeoData().GetEllipseAxes()));
                                DivideOverAxis(cell, axis);
                        }
                        division_count++;
                }
        }
        return division_count;
}

void CellDivider::DivideOverAxis(Cell* cell, const array<double, 3>& axis)
{
        assert((m_cd.m_mesh != nullptr) && "Mesh not initialized.");

        //--------------------------------------------------------------------------
        // First we need to find which edges get intersected by the division axis.
        // The algorithm detects a change of sign in the cross product of the
        // axis and a vector form to centroid to the node. If the sign flips, the
        // edge from the previous node to the current node gets intersected. We
        // store the iterators to the current node in the intersect_nodes.
        //--------------------------------------------------------------------------
        auto intersect_nodes = ComputeIntersectNodes(cell, axis);
        assert((intersect_nodes.size() == 2) && "Can only handle two intersection points.");

        //--------------------------------------------------------------------------
        // Division currently only works for two intersection points only. When there
        // are more than two intersection points, the two points that are closest to the
        // centroid on either side of the long axis of the cell are chosen as the
        // intersection points.
        // The division creates a new cell (one daughter) and morphs the parent
        // into a new cell (the other daughter). We need to manage the parent's
        // attributes (cell and some geometric properties) to be able to assign
        // attributes for each of the daughters.
        //--------------------------------------------------------------------------
        if (intersect_nodes.size() != 2) {
                // ----------------------------------------------------------------
                // If there are less than two intersection points or the centroid
                // is not in the interior of the cell, still throw exception
                // ----------------------------------------------------------------
                if((Geom::wn_PnPoly(cell->GetCentroid(), cell->GetNodes()) != 0)
                        && intersect_nodes.size() > 2)
                {
                        intersect_nodes = ComputeTwoIntersectNodes(cell, axis, intersect_nodes);
                } else {
                        const string msg = string(VL_HERE) + " exception:\nToo many intersection points.";
                        throw Exception(msg);
                }
        }

        //--------------------------------------------------------------------------
        // Save parent's attributes to fix the daughter's after split.
        //--------------------------------------------------------------------------
        CellAttributes  this_cell_old_attributes = *cell;

        //--------------------------------------------------------------------------
        // Split geometry of parent (cell prior to split) along a line through the
        // centroid in the direction of the axis to obtain the daughters (resp. cell
        // after split and cell in return value).
        //--------------------------------------------------------------------------
        const auto from = cell->GetCentroid();
        const auto to   = from + axis;
        Cell* daughter = SplitCell(cell, intersect_nodes, from, to);
        cell->IncrementDivisionCounter();
        daughter->IncrementDivisionCounter();

        //--------------------------------------------------------------------------
        // Update daughter attributes.
        //--------------------------------------------------------------------------
        daughter->CellAttributes::operator=(this_cell_old_attributes);
        // solute
        const auto new_solute       = cell->GetSolute() / sqrt(2.0);
        cell->SetSolute(new_solute);
        daughter->SetSolute(new_solute);
        // target area
        const auto new_t_area       = cell->GetTargetArea() / 2.0;
        cell->SetTargetArea(new_t_area);
        daughter->SetTargetArea(new_t_area);
        // target length
        const auto new_t_length     = cell->GetTargetLength() / sqrt(2.0);
        cell->SetTargetLength(new_t_length);
        daughter->SetTargetLength(new_t_length);
        // Chemicals in cells and transporters in walls.
        m_cell_daughters(daughter, cell);
}

list<Cell*> CellDivider::GeometricSplitCell(Cell* cell, Node* node1, Node* node2)
{
        vector<Node*> intersect_nodes;
        vector<Node*> nodes = cell->GetNodes();
        for (unsigned int i = 0; i < nodes.size(); i++) {
                if (nodes.at(i) == node1 || nodes.at(i) == node2) {
                        intersect_nodes.push_back(nodes.at(i));
                }
        }

        Cell* newCell = SplitCell(cell, intersect_nodes, *node1, *node2);
        cell->SetGeoDirty();
        newCell->SetGeoDirty();
        newCell->CellAttributes::operator=(*cell);

        list<Cell*> newCells;
        newCells.push_back(cell);
        newCells.push_back(newCell);

        return newCells;
}

void CellDivider::Initialize(Mesh* mesh)
{
        try {
                assert( (mesh != nullptr) && "mesh pointer not ok.");
                m_cd                       = CoreData();
                m_collapse_node_threshold  = 10e-12;    //Due to floating point errors.
                m_copy_wall                = true;
                m_mesh                     = mesh;
                m_model_name               = "";
                m_target_node_distance     = 0.0;

        }
        catch(exception& e) {
                const string here = string(VL_HERE) + " exception:\n";
                throw Exception(here + e.what());
        }
}

void CellDivider::Initialize(const CoreData& cd)
{
        try {
                assert( cd.Check() && "CoreData not ok.");
                m_cd                       = cd;
                auto& p                    = m_cd.m_parameters;
                m_mesh                     = m_cd.m_mesh.get();
                const auto model_group     = m_cd.m_parameters->get<string>("model.group", "");

                const auto factory = ComponentFactoryProxy::Create(model_group);
                m_cell_daughters = factory->CreateCellDaughters(m_cd);
                m_cell_split = factory->CreateCellSplit(m_cd);

                m_collapse_node_threshold  = 0.05;
                m_copy_wall                = p->get<bool>("cell_mechanics.copy_wall");
                m_model_name               = p->get<string>("model.name");
                m_target_node_distance     = p->get<double>("cell_mechanics.target_node_distance");

                const trng::uniform01_dist<double> dist;
                m_uniform_generator        = m_cd.m_random_engine->GetGenerator(dist);
        }
        catch(exception& e) {
                const string here = string(VL_HERE) + " exception:\n";
                throw Exception(here + e.what());
        }
}

Node* CellDivider::InsertNewNodeAtIntersect(Cell* this_cell, Cell* neighbor_cell,
        array<double, 3> node_vec, const Edge& intrsct_edge)
{
        auto& this_cell_nodes   = this_cell->GetNodes();
        const bool at_boundary  = m_mesh->IsAtBoundary(&intrsct_edge);
        Node* node_ptr          = m_mesh->BuildNode(node_vec, at_boundary);

        //--------------------------------------------------------------------------
        // Insert new node in this cell in the appropriate position.
        //--------------------------------------------------------------------------
        const auto ins_pos_this_cell = find(this_cell_nodes.begin(), this_cell_nodes.end(), intrsct_edge.GetSecond());
        assert(((ins_pos_this_cell != this_cell_nodes.begin() && *prev(ins_pos_this_cell) == intrsct_edge.GetFirst())
                || this_cell_nodes.back() == intrsct_edge.GetFirst()) && "Error in insertion position!");
        this_cell_nodes.insert(ins_pos_this_cell, node_ptr);

        //--------------------------------------------------------------------------
        // Insert new node in neighbor cell. First, find the position of the first node.
        //--------------------------------------------------------------------------
        auto& nb_nodes = neighbor_cell->GetNodes();
#if (defined(__clang__) && defined(__APPLE__))
        const auto ins_pos_nb_cell = find(nb_nodes.cbegin(), nb_nodes.cend(), intrsct_edge.GetFirst());
#else
        const auto ins_pos_nb_cell = find(nb_nodes.begin(), nb_nodes.end(), intrsct_edge.GetFirst());
#endif

#if (defined(__clang__) && defined(__APPLE__))
        const auto cit = make_const_circular(nb_nodes.cbegin(), nb_nodes.cend(), ins_pos_nb_cell);
#else
        const auto cit = make_const_circular(nb_nodes.begin(), nb_nodes.end(), ins_pos_nb_cell);
#endif
        //--------------------------------------------------------------------------
        // The node comes before or after the first (depending on whether the neighbor
        // cell is a real cell or the boundary polygon this differs: you cannot simply
        // deduce this by reversing the order from what it was in ''this cell''.
        //--------------------------------------------------------------------------
        if (*next(cit) == intrsct_edge.GetSecond()) {
                nb_nodes.insert((next(cit)).get(), node_ptr);
        } else {
                assert((*prev(cit) == intrsct_edge.GetSecond()) && "Error in insertion position!");
                nb_nodes.insert(cit.get(), node_ptr);
        }

        return node_ptr;
}

Cell* CellDivider::SplitCell(Cell* this_cell, vector<Node*> intersect_nodes,
        const array<double, 3>& from, const array<double, 3>& to)
{
        auto daughter               = m_mesh->BuildCell();
        auto& this_cell_nodes       = this_cell->GetNodes();
        auto& daughter_nodes        = daughter->GetNodes();

        array<array<double, 3>, 2>  new_node_vector;
        array<Node*, 2>             new_node_ptr;
        array<int, 2>               new_node_flag;
        array<Edge, 2>              intersect_edge;
        array<Wall*, 4>             div_wall    {{nullptr, nullptr, nullptr, nullptr}};
        array<double, 2>            orig_length {{0.0, 0.0}};
        array<double, 2>            orig_rest_length {{0.0, 0.0}};
        array<double, 2>            orig_rest_length_init {{0.0, 0.0}};

        //--------------------------------------------------------------------------
        // Preliminary: compute required info to effect the split: which edges are
        // intersected, what are the intersection points, do the intersection
        // points need to be coalesced with existing nodes because the are so
        // close to them.
        //--------------------------------------------------------------------------
        ComputeIntersection(this_cell, intersect_nodes, from, to,
                new_node_vector, new_node_ptr, new_node_flag, intersect_edge);

        //--------------------------------------------------------------------------
        // For both edges that have an intersection on them: insert its new node
        // into all cells that own the edge but only if it really is a new node.
        // If the intersect is too close to an existing node, do nothing.
        //--------------------------------------------------------------------------
        for (int i = 0; i < 2; i++) {
                Cell* neighbor_cell = m_mesh->FindEdgeNeighbor(this_cell, intersect_edge[i]);
                assert(neighbor_cell != nullptr && "Where is the neighbor?");

                if (new_node_flag[i] == 0) {
                        new_node_ptr[i] = InsertNewNodeAtIntersect(
                                this_cell, neighbor_cell, new_node_vector[i], intersect_edge[i]);

                        //--------------------------------------------------------------------------
                        // Update node ownership.
                        //--------------------------------------------------------------------------
                        m_mesh->UpdateNodeOwningNeighbors(neighbor_cell);

                        //--------------------------------------------------------------------------
                        // If node is inserted into fixed edge (i.e. in the petiole)
                        // make the new node fixed as well.
                        //--------------------------------------------------------------------------
                        (new_node_ptr[i])->SetFixed(intersect_edge[i].IsFixed());
                }

                AddSplitWalls(this_cell, neighbor_cell, i, new_node_ptr, div_wall, orig_length, orig_rest_length, orig_rest_length_init);

        }

        //--------------------------------------------------------------------------
        // Bookkeeping for the nodes:
        // 1) existing nodes need to be re-assigned from parent to the daughters.
        // 2) add nodes of what will be the shared wall i.e those of the intersection
        //    plus intermediate ones along the wall if necessary
        // 3) repair cell <-> node connectivities
        //--------------------------------------------------------------------------
        AssignExistingNodesToSplitCells(
                m_mesh, this_cell, this_cell_nodes, daughter_nodes, new_node_ptr);
        AddSharedWallNodes(
                new_node_vector, new_node_ptr, this_cell_nodes, daughter_nodes);
        m_mesh->UpdateNodeOwningNeighbors(this_cell);
        m_mesh->UpdateNodeOwningNeighbors(daughter);

        //--------------------------------------------------------------------------
        // Recalculate inertia integrals (this includes area)
        //--------------------------------------------------------------------------
        this_cell->SetGeoDirty();
        daughter->SetGeoDirty();

        //--------------------------------------------------------------------------
        // Bookkeeping for the walls:
        // Existing walls need to reassigned and the cell <-> wall connectivity updated.
        // The new wall between daughters is not included (need not be transferred).
        //--------------------------------------------------------------------------
        auto copy_walls = this_cell->GetWalls();
        for (Wall* w : copy_walls) {
                auto first = find(this_cell_nodes.begin(), this_cell_nodes.end(), w->GetN1());
                auto second = find(this_cell_nodes.begin(), this_cell_nodes.end(), w->GetN2());

                if (first == this_cell_nodes.end() || second == this_cell_nodes.end()
                        || (w->GetN1() == new_node_ptr[0] && w->GetN2() == new_node_ptr[1])) {

                        this_cell->ReassignWall(w, daughter);
                        if (w->GetC1() == this_cell) {
                                w->SetC1(daughter);
                        } else {
                                w->SetC2(daughter);
                        }
                }
        }

        //--------------------------------------------------------------------------
        // Special case: a single cell is split. Single cells have no walls at all
        // so we need to build the outer walls first.
        //--------------------------------------------------------------------------
        if (m_mesh->GetCells().size() == 2) {
                Wall* daughter_wall = m_mesh->BuildWall(
                        new_node_ptr[0], new_node_ptr[1], daughter,
                        m_mesh->GetBoundaryPolygon());
                daughter->AddWall(daughter_wall);
                m_mesh->GetBoundaryPolygon()->AddWall(daughter_wall);
                // Correct the walls associated with the nodes on the wall
                m_mesh->UpdateNodeOwningWalls(daughter_wall);

                Wall* parent_wall = m_mesh->BuildWall(
                        new_node_ptr[1], new_node_ptr[0], this_cell,
                        m_mesh->GetBoundaryPolygon());
                this_cell->AddWall(parent_wall);
                m_mesh->GetBoundaryPolygon()->AddWall(parent_wall);
                // Correct the walls associated with the nodes on the wall
                m_mesh->UpdateNodeOwningWalls(parent_wall);
        }

        //--------------------------------------------------------------------------
        // Add the new wall shared by the daughter and to the necessary bookkeeping.
        //--------------------------------------------------------------------------
        auto wall = m_mesh->BuildWall(
                new_node_ptr[0], new_node_ptr[1], this_cell, daughter);
        this_cell->AddWall(wall);
        daughter->AddWall(wall);
        // Correct the walls associated with the nodes on the wall
        m_mesh->UpdateNodeOwningWalls(wall);

        if ( m_model_name == "Wortel" ) {
                double wall_strength = 0.0;
                std::array<double, 3> ref{{1., 0., 0.}};
                std::array<double, 3> tot1 = *wall->GetN1();
                std::array<double, 3> tot2 = *wall->GetN2();
                std::array<double, 3> tot3 = tot1 - tot2;
                double t2 = SignedAngle(tot3, ref);
                if ( ( t2 > (-pi() * 3 / 4) && t2 <= (-pi() / 4) ) || ( t2 > (pi() / 4) && t2 <= (pi() * 3 / 4) ) )
                {
                        wall_strength = 1.;
                } else {
                        wall_strength = 5.;//DDV: should be 1 for Blad, 5 for Wortel.
                }
                wall->SetStrength(wall_strength);
        }

        //--------------------------------------------------------------------------
        // Correct transporter concentrations of divided walls.
        //--------------------------------------------------------------------------
        for (int i = 0; i < 4; i++) {
                if (div_wall[i]) {
                        const double current_length = div_wall[i]->GetLength();
                        const double correction = current_length / orig_length[i/2];

                        div_wall[i]->SetRestLength(orig_rest_length[i/2] * correction);
                        div_wall[i]->SetRestLengthInit(orig_rest_length_init[i/2] * correction);

                        auto transporters1 = div_wall[i]->GetTransporters1();
                        for (auto& t : transporters1) {
                                t *= correction;
                        }
                        div_wall[i]->SetTransporters1(transporters1);

                        auto transporters2 = div_wall[i]->GetTransporters2();
                        for (auto& t : transporters2) {
                                t *= correction;
                        }
                        div_wall[i]->SetTransporters2(transporters2);
                }
        }

        //--------------------------------------------------------------------------
        // Collect neighbors affected by this cell's division, namely
        // this cell's old neighbors, this cell and the daughter.
        //--------------------------------------------------------------------------
        list<Cell*>  broken_neighbors;
        auto&        this_neighbors = m_mesh->GetNeighbors(this_cell);
        copy(this_neighbors.begin(), this_neighbors.end(), back_inserter(broken_neighbors));
        broken_neighbors.push_back(this_cell);
        broken_neighbors.push_back(daughter);

        //--------------------------------------------------------------------------
        // Now reconstruct neighbor list for all "broken" neighbors.
        //--------------------------------------------------------------------------
        for (auto const& bn : broken_neighbors) {
                m_mesh->ConstructNeighborList(bn);
        }

        return daughter;
}

} // namespace