EditableMesh.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 "EditableMesh.h"

#include "algo/CellDivider.h"
#include "bio/Cell.h"
#include "bio/CellAttributes.h"
#include "bio/Edge.h"
#include "bio/Mesh.h"
#include "bio/MeshCheck.h"
#include "bio/NeighborNodes.h"
#include "bio/Node.h"
#include "bio/PBMBuilder.h"
#include "bio/Wall.h"
#include "bio/WallAttributes.h"
#include "editor/PolygonUtils.h"
#include "util/container/circular_iterator.h"

#include <algorithm>

namespace SimPT_Sim { class Mesh; }

namespace SimPT_Editor {

using namespace std;
using namespace SimPT_Sim;
using namespace SimPT_Sim::Container;

EditableMesh::EditableMesh(const boost::property_tree::ptree& pt)
        : m_mesh(PBMBuilder().Build(pt)), m_cells(m_mesh->GetCells()),
          m_nodes(m_mesh->GetNodes()), m_walls(m_mesh->GetWalls()) {}

EditableMesh::~EditableMesh() {}

bool EditableMesh::CanSplitCell(Node* node1, Node* node2, Cell* cell) const
{
        list<Cell*> owners;
        for (auto nb : m_mesh->GetNodeOwningNeighbors(node1)) {
                if (nb.GetCell()->GetIndex() != -1) {
                        owners.push_back(nb.GetCell());
                }
        }
        for (auto nb : m_mesh->GetNodeOwningNeighbors(node2)) {
                if (nb.GetCell()->GetIndex() != -1) {
                        owners.push_back(nb.GetCell());
                }
        }
        owners.sort([](Cell* c1, Cell* c2) {return c1 < c2;});

        auto eql = [](Cell* c1, Cell* c2) {
                return c1 == c2;
        };

        list<Cell*> commonOwners;
        auto owner = adjacent_find(owners.begin(), owners.end(), eql);
        while (owner != owners.end()) {
                commonOwners.push_back(*owner);
                owner = adjacent_find(++owner, owners.end(), eql);
        }

        if (commonOwners.empty()) {
                return false;
        }
        else if (cell != nullptr) {
                if (find(commonOwners.begin(), commonOwners.end(), cell) != commonOwners.end()) {
                        commonOwners.clear();
                        commonOwners.push_back(cell);
                }
                else {
                        return false;
                }
        }

        QLineF cut((*node1)[0], (*node1)[1], (*node2)[0], (*node2)[1]);
        for (Cell* c : commonOwners) {
                QPolygonF polygon;
                for (Node* n : c->GetNodes()) {
                        polygon.push_back(QPointF((*n)[0], (*n)[1]));
                }

                auto polygons = PolygonUtils::SlicePolygon(polygon, cut);
                if (polygons.size() == 2) {
                        bool ok = true;
                        auto it = polygons.begin();
                        while (it != polygons.end() && ok) {
                                if (!it->contains(cut.p1()) || !it->contains(cut.p2())) {
                                        ok = false;
                                }
                                for (const QPointF& p : *it) {
                                        if (!polygon.contains(p)) {
                                                ok = false;
                                        }
                                }
                                it++;
                        }
                        if (ok) {
                                return true;
                        }
                }
        }

        return false;
}

Cell* EditableMesh::CreateCell(Node* node1, Node* node2, const QPointF& newPoint)
{
        assert(node1->IsAtBoundary() && node2->IsAtBoundary()
                && "The two given nodes aren't located at the boundary of the mesh.");
        assert(find(m_mesh->GetBoundaryPolygon()->GetNodes().begin(), m_mesh->GetBoundaryPolygon()->GetNodes().end(), node1) != m_mesh->GetBoundaryPolygon()->GetNodes().end()
                && find(m_mesh->GetBoundaryPolygon()->GetNodes().begin(), m_mesh->GetBoundaryPolygon()->GetNodes().end(), node1) != m_mesh->GetBoundaryPolygon()->GetNodes().end()
                && "The two given nodes aren't located at the boundary of the mesh.");

        QPolygonF boundaryPoints;
        for (Node* n : m_mesh->GetBoundaryPolygon()->GetNodes()) {
                boundaryPoints.append(QPointF((*n)[0], (*n)[1]));
        }

        list<QLineF> excluded;
        list<QLineF> included({QLineF(QPointF((*node1)[0], (*node1)[1]), newPoint), QLineF(QPointF((*node2)[0], (*node2)[1]), newPoint)});
        if (boundaryPoints.contains(newPoint) || boundaryPoints.containsPoint(newPoint, Qt::OddEvenFill)
                || !IsConsistent(excluded, included)) {

                return nullptr;
        }
        else {
                // There are two possible polygons which can be the cell.
                // The two ways you can go around the boundary polygon.
                list<QPolygonF> polygons({QPolygonF(), QPolygonF()});
                auto current = polygons.begin();

                QPointF point1((*node1)[0], (*node1)[1]);
                QPointF point2((*node2)[0], (*node2)[1]);

                for (const QPointF& p : boundaryPoints) {
                        current->append(p);
                        if (p == point1 || p == point2) {
                                current->append(newPoint);
                                if (++current == polygons.end()) {
                                        current = polygons.begin();
                                }
                                current->append(p);
                        }
                }

                // The smallest one is the right one, since the largest one will form the boundary polygon.
                QPolygonF polygon = polygons.front().intersected(polygons.back());
                if (polygon.front() == polygon.back()) {
                        polygon.pop_back(); //So the polygon will become open.
                }

                // Make the polygon counterclockwise (in the Qt coordinate system).
                if (PolygonUtils::IsClockwise(polygon)) {
                        reverse(polygon.begin(), polygon.end());
                }

                // After this, the nodes are in the order in which they will appear in the wall
                // between the new cell and the boundary polygon.
                rotate(polygon.begin(), find(polygon.begin(), polygon.end(), newPoint), polygon.end());
                if (*next(polygon.begin()) == point1) {
                        swap(node1, node2);
                        swap(point1, point2);
                }
                else {
                        assert(*next(polygon.begin()) == point2
                                && "One of the two points should be next to the new point.");
                }

                //Build the new node.
                auto newNode = m_mesh->BuildNode({{newPoint.x(), newPoint.y(), 0}}, true);

                //Order old boundary nodes to have the first node (start of the wall) at the front.
                auto oldBoundaryNodes = m_mesh->GetBoundaryPolygon()->GetNodes();
                auto node_it = find(oldBoundaryNodes.begin(), oldBoundaryNodes.end(), node1);
                rotate(oldBoundaryNodes.begin(), node_it, oldBoundaryNodes.end());

                //Create the new boundary nodes with the first node (start of the wall) at the front.
                auto& newBoundaryNodes = m_mesh->GetBoundaryPolygon()->GetNodes();
                node_it = find(newBoundaryNodes.begin(), newBoundaryNodes.end(), node2);
                rotate(newBoundaryNodes.begin(), node_it, newBoundaryNodes.end());
                node_it = find(newBoundaryNodes.begin(), newBoundaryNodes.end(), node1);
                newBoundaryNodes.erase(next(node_it), newBoundaryNodes.end());
                newBoundaryNodes.push_back(newNode); // may invalidate node_it, so compute it again
                node_it = find(newBoundaryNodes.begin(), newBoundaryNodes.end(), node1);
                rotate(newBoundaryNodes.begin(), node_it, newBoundaryNodes.end());
                m_mesh->GetBoundaryPolygon()->SetGeoDirty();

                // Update the walls.
                if (m_walls.size() == 0) {
                        assert(m_cells.size() == 1 && "There should be two cells in the mesh.");

                        Cell* cell = m_cells.front();

                        Wall* newWall = m_mesh->BuildWall(node2, node1, cell, m_mesh->GetBoundaryPolygon());
                        newWall->GetC1()->AddWall(newWall);
                        newWall->GetC2()->AddWall(newWall);
                        m_mesh->UpdateNodeOwningWalls(newWall);

                        newWall = m_mesh->BuildWall(node1, node2, cell, m_mesh->GetBoundaryPolygon());
                        newWall->GetC1()->AddWall(newWall);
                        newWall->GetC2()->AddWall(newWall);
                        m_mesh->UpdateNodeOwningWalls(newWall);
                }
                else {
                        // This case requires the possible splitting of walls.
                        assert(m_walls.size() >= 2 && "There can't be one wall in the mesh.");

                        auto& boundaryWalls = m_mesh->GetBoundaryPolygon()->GetWalls();

                        if (m_mesh->GetNodeOwningNeighbors(node1).size() == 2) {
                                auto endpoint = oldBoundaryNodes.rbegin();
                                while (m_mesh->GetNodeOwningNeighbors(*endpoint).size() == 2) {
                                        endpoint++;
                                        assert(endpoint != oldBoundaryNodes.rend()
                                                && "There should be at least two endpoints with more than 2 owners in the mesh.");
                                }

                                auto wall = find_if(boundaryWalls.begin(), boundaryWalls.end(), [endpoint](Wall* w){return w->GetN1() == *endpoint;});
                                assert(wall != boundaryWalls.end() && "There should be a wall for this endpoint.");

                                Wall* newWall = m_mesh->BuildWall(node1, (*wall)->GetN2(), (*wall)->GetC1(), (*wall)->GetC2());
                                newWall->WallAttributes::operator=(**wall);
                                newWall->GetC1()->AddWall(newWall);
                                newWall->GetC2()->AddWall(newWall);

                                (*wall)->SetN2(node1);
                                (*wall)->SetLengthDirty();

                                m_mesh->UpdateNodeOwningWalls(newWall);
                        }

                        if (m_mesh->GetNodeOwningNeighbors(node2).size() == 2) {
                                auto endpoint = find(oldBoundaryNodes.begin(), oldBoundaryNodes.end(), node2);
                                while (m_mesh->GetNodeOwningNeighbors(*endpoint).size() == 2) {
                                        endpoint++;
                                        assert(endpoint != oldBoundaryNodes.end() && "There should be at least two endpoints with more than 2 owners in the mesh.");
                                }

                                auto wall = find_if(boundaryWalls.begin(), boundaryWalls.end(), [endpoint](Wall* w){return w->GetN2() == *endpoint;});
                                assert(wall != boundaryWalls.end() && "There should be a wall for this endpoint.");

                                Wall* newWall = m_mesh->BuildWall((*wall)->GetN1(), node2, (*wall)->GetC1(), (*wall)->GetC2());
                                newWall->WallAttributes::operator=(**wall);
                                newWall->GetC1()->AddWall(newWall);
                                newWall->GetC2()->AddWall(newWall);

                                (*wall)->SetN1(node2);
                                (*wall)->SetLengthDirty();

                                m_mesh->UpdateNodeOwningWalls(newWall);
                        }
                }

                // Make a vector of node IDs for the new cell (in the right order).
                list<Node*> cellNodes(oldBoundaryNodes.begin(), next(find(oldBoundaryNodes.begin(), oldBoundaryNodes.end(), node2)));
                cellNodes.reverse();
                cellNodes.push_front(newNode);

                vector<unsigned int> nodeIDs(cellNodes.size());
                transform(cellNodes.begin(), cellNodes.end(), nodeIDs.begin(), [](Node* n){return n->GetIndex();});

                // Build the new cell.
                auto newCell = m_mesh->BuildCell(nodeIDs);

                // Reassign the correct walls from the boundary polygon to the new cell.
                const auto boundaryWalls = m_mesh->GetBoundaryPolygon()->GetWalls();
                for (Wall* w : boundaryWalls) {
                        auto n1 = find(cellNodes.begin(), cellNodes.end(), w->GetN1());
                        auto n2 = find(cellNodes.begin(), cellNodes.end(), w->GetN2());
                        if (n1 != cellNodes.end() && n2 != cellNodes.end() && (*n2 != node1 || *n1 != node2)) {
                                w->SetC2(newCell);
                                m_mesh->GetBoundaryPolygon()->ReassignWall(w, newCell);
                        }
                }

                // Create the new wall between the new cell and the boundary polygon.
                Wall* newWall = m_mesh->BuildWall(node1, node2,
                        newCell, m_mesh->GetBoundaryPolygon());
                newWall->GetC1()->AddWall(newWall);
                newWall->GetC2()->AddWall(newWall);
                m_mesh->UpdateNodeOwningWalls(newWall);

                // Update the boundary neighbors of the nodes in the new cell to have the correct neighbors.
                for (Node* n : newCell->GetNodes()) {
                        if (n != newNode) {
                                auto& owners = m_mesh->GetNodeOwningNeighbors(n);
                                auto owner = find_if(owners.begin(), owners.end(),
                                        [](const NeighborNodes& nb){return nb.GetCell()->GetIndex() == -1;});
                                assert(owner != owners.end()
                                        && "Node without a boundary neighbor in new cell that wasn't new node.");
                                owners.erase(owner);
                        }
                }

                auto boundary_cit = make_const_circular(newBoundaryNodes);
                do {
                        m_mesh->GetNodeOwningNeighbors(*boundary_cit).push_back(
                                NeighborNodes(m_mesh->GetBoundaryPolygon(),
                                                *prev(boundary_cit), *next(boundary_cit)));
                } while (*boundary_cit++ != node2);

                //Construct the neighbor lists of the new cell.
                m_mesh->ConstructNeighborList(newCell);
                for (Cell* c : m_mesh->GetNeighbors(newCell)) {
                        m_mesh->ConstructNeighborList(c);
                }
                m_mesh->ConstructNeighborList(m_mesh->GetBoundaryPolygon());

                // Set the nodes of the cell to the boundary or not when necessary.
                for (Node* n : newCell->GetNodes()) {
                        SetAtBoundary(n);
                }

                assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");
                return newCell;
        }
}

void EditableMesh::DeleteCells(const list<Cell*>& cells)
{
        assert(cells.size() < m_cells.size() && "Can't delete all cells in the mesh.");
        assert(find_if(cells.begin(), cells.end(), [this](Cell* c) {return (find(m_cells.begin(), m_cells.end(), c) == m_cells.end());}) == cells.end()
                && "There's a given cell that doesn't exist in the mesh.");

        if (cells.empty()) {
                return;
        }

        // Delete the given cells and update the walls.
        auto cell = m_cells.begin();
        while (cell != m_cells.end()) {
                if (find(cells.begin(), cells.end(), *cell) != cells.end()) {
                        for (Node* n : (*cell)->GetNodes()) {
                                n->SetAtBoundary(true);
                        }

                        for (Wall* w : (*cell)->GetWalls()) {
                                if (w->IsAtBoundary()) {
                                        w->SetN1(nullptr);
                                        w->SetN2(nullptr);
                                        m_walls.erase(remove_if(m_walls.begin(), m_walls.end(),
                                                        [w](Wall* wall){return wall == w;}),
                                                        m_walls.end());
                                        //m_walls.remove_if([w](Wall* wall){return wall == w;});
                                }
                                else if (w->GetC1() == *cell) {
                                        w->SetC1(w->GetC2());
                                        w->SetC2(m_mesh->GetBoundaryPolygon());

                                        Node* tmpNode = w->GetN1();
                                        w->SetN1(w->GetN2());
                                        w->SetN2(tmpNode);

                                        // No need to set the wall length dirty bit, as this is the same wall:
                                        // The cells and nodes are just swapped so that C1 isn't the boundary polygon
                                }
                                else { //walls->GetC2() == cell
                                        w->SetC2(m_mesh->GetBoundaryPolygon());
                                }
                        }

                        m_mesh->GetCells().erase(cell);
                }
                else {
                        cell++;
                }
        }

        //If there is only one cell left, all walls need to be deleted.
        if (m_cells.size() == 1) {
                m_walls.clear();
                m_cells.front()->GetWalls().clear();
                m_mesh->GetBoundaryPolygon()->GetWalls().clear();
        }

        for (auto n : m_nodes) {
                m_mesh->GetNodeOwningNeighbors(n).clear();
        }

        list<Node*> nodesToKeep;
        for (auto c : m_cells) {
                auto cit = make_const_circular(c->GetNodes(), c->GetNodes().begin());
                do {
                        nodesToKeep.push_back(*cit);
                        m_mesh->GetNodeOwningNeighbors(*cit).push_back(
                                NeighborNodes(c, *prev(cit), *next(cit)));
                } while (++cit != c->GetNodes().begin());
        }

        // Delete the useless nodes.
        list<Node*> boundaryNodes;
        auto node = m_nodes.begin();
        while (node != m_nodes.end()) {
                if (find(nodesToKeep.begin(), nodesToKeep.end(), *node) == nodesToKeep.end()) {
                        m_mesh->RemoveNodeFromOwnerLists(*node);
                        m_nodes.erase(node);
                }
                else {
                        if ((*node)->IsAtBoundary()) {
                                boundaryNodes.push_back(*node);
                        }
                        node++;
                }
        }

        //Re-add the new boundaryNodes to the boundary polygon.
        Node* current = boundaryNodes.front();
        boundaryNodes.pop_front();
        m_mesh->GetBoundaryPolygon()->GetNodes().clear();
        m_mesh->GetBoundaryPolygon()->GetNodes().push_back(current);
        while (!boundaryNodes.empty()) {
                // These are to check whether there are two neighboring
                // cells (which implies we have the wrong boundary node).
                list<Node*> firstNeighbors;
                // These are the real options, as the order of the nodes
                // should be the same for all the cells (counterclockwise).
                list<Node*> secondNeighbors;
                for (NeighborNodes nb : m_mesh->GetNodeOwningNeighbors(current)) {
                        if (nb.GetNb1()->IsAtBoundary()) {
                                firstNeighbors.push_back(nb.GetNb1());
                        }
                        if (nb.GetNb2()->IsAtBoundary()) {
                                secondNeighbors.push_back(nb.GetNb2());
                        }
                }

                auto next = find_if(boundaryNodes.begin(), boundaryNodes.end(),
                        [firstNeighbors, secondNeighbors](Node* n) {
                                return find(firstNeighbors.begin(), firstNeighbors.end(), n)
                                                == firstNeighbors.end()
                                        && find(secondNeighbors.begin(), secondNeighbors.end(), n)
                                                != secondNeighbors.end();
                });
                assert(next != boundaryNodes.end()
                        && "Inconsistency: No boundary neighbor to current boundary node.");

                current = *next;
                boundaryNodes.erase(next);
                m_mesh->GetBoundaryPolygon()->GetNodes().push_back(current);
        }

        auto cit = make_const_circular(m_mesh->GetBoundaryPolygon()->GetNodes(),
                                m_mesh->GetBoundaryPolygon()->GetNodes().begin());
        do {
                m_mesh->GetNodeOwningNeighbors(*cit).push_back(
                        NeighborNodes(m_mesh->GetBoundaryPolygon(), *prev(cit), *next(cit)));
        } while (++cit != m_mesh->GetBoundaryPolygon()->GetNodes().begin());

        // Boundary walls.
        auto cmp_boundary_walls = [](Wall* w1, Wall* w2) {
                Cell* c1 = w1->GetC1()->GetIndex() != -1 ? w1->GetC1() : w1->GetC2();
                Cell* c2 = w2->GetC1()->GetIndex() != -1 ? w2->GetC1() : w2->GetC2();
                return c1->GetIndex() < c2->GetIndex();
        };
        auto eql_boundary_walls = [](Wall* w1, Wall* w2) {
                Cell* c1 = w1->GetC1()->GetIndex() != -1 ? w1->GetC1() : w1->GetC2();
                Cell* c2 = w2->GetC1()->GetIndex() != -1 ? w2->GetC1() : w2->GetC2();
                return c1->GetIndex() == c2->GetIndex();
        };

        list<Wall*> boundaryWalls;
        for (auto w : m_walls) {
                if (w->IsAtBoundary()) {
                        boundaryWalls.push_back(w);
                }
        }

        //TODO: Random merge of walls. Something better?
        m_mesh->GetBoundaryPolygon()->GetWalls().clear();
        boundaryWalls.sort(cmp_boundary_walls);
        auto wall = adjacent_find(boundaryWalls.begin(), boundaryWalls.end(), eql_boundary_walls);
        while (wall != boundaryWalls.end()) {
                auto next = wall;
                next++;
                while (next != boundaryWalls.end() && eql_boundary_walls(*wall, *next)) {
                        bool merge = true;
                        if ((*wall)->GetN1() == (*next)->GetN1()) {
                                (*wall)->SetN1((*next)->GetN2());
                        }
                        else if ((*wall)->GetN1() == (*next)->GetN2()) {
                                (*wall)->SetN1((*next)->GetN1());
                        }
                        else if ((*wall)->GetN2() == (*next)->GetN1()) {
                                (*wall)->SetN2((*next)->GetN2());
                        }
                        else if ((*wall)->GetN2() == (*next)->GetN2()) {
                                (*wall)->SetN2((*next)->GetN1());
                        }
                        else {
                                merge = false;
                        }

                        if (merge) {
                                (*wall)->SetLengthDirty();

                                (*next)->GetC1()->GetWalls().remove(*next);
                                (*next)->GetC2()->GetWalls().remove(*next);

                                Wall* tmpWall = *next;
                                boundaryWalls.erase(next);

                                //Remove from the node owning walls list.
                                tmpWall->SetN1(nullptr);
                                tmpWall->SetN2(nullptr);
                                m_mesh->UpdateNodeOwningWalls(*wall);

                                //m_walls.remove_if([tmpWall](Wall* w){return w == tmpWall;});
                                m_walls.erase(remove_if(
                                        m_walls.begin(), m_walls.end(),
                                        [tmpWall](Wall* w){return w == tmpWall;}),
                                        m_walls.end());
                                next = wall;
                        }
                        next++;
                }
                wall = adjacent_find(++wall, boundaryWalls.end(), eql_boundary_walls);
        }

        for (auto w : boundaryWalls) {
                m_mesh->GetBoundaryPolygon()->GetWalls().push_back(w);
        }

        m_mesh->GetBoundaryPolygon()->SetGeoDirty();

        m_mesh->ConstructNeighborList(m_mesh->GetBoundaryPolygon());

        for (auto c : m_mesh->GetCells()) {
                m_mesh->ConstructNeighborList(c);
        }

        // Fix the IDs.
        FixNodeIDs();
        FixWallIDs();
        FixCellIDs();

        if (m_mesh->GetWalls().size() == 0) {
                for (auto n : m_mesh->GetNodes()) {
                        m_mesh->GetNodeOwningWalls(n).clear();
                }
        }

        for (auto w : m_mesh->GetWalls()) {
                m_mesh->UpdateNodeOwningWalls(w);
        }

        assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");
}

bool EditableMesh::DeleteTwoDegreeNode(Node* node)
{
        list<QLineF> oldEdges;
        list<QLineF> newEdges;
        for (NeighborNodes neighbor : m_mesh->GetNodeOwningNeighbors(node)) {
                QLineF oldEdge1((*neighbor.GetNb1())[0], (*neighbor.GetNb1())[1], (*node)[0], (*node)[1]);
                QLineF oldEdge2((*neighbor.GetNb2())[0], (*neighbor.GetNb2())[1], (*node)[0], (*node)[1]);
                QLineF newEdge((*neighbor.GetNb1())[0], (*neighbor.GetNb1())[1], (*neighbor.GetNb2())[0], (*neighbor.GetNb2())[1]);

                if (find(oldEdges.begin(), oldEdges.end(), oldEdge1) == oldEdges.end()) {
                        oldEdges.push_front(oldEdge1);
                }
                if (find(oldEdges.begin(), oldEdges.end(), oldEdge2) == oldEdges.end()) {
                        oldEdges.push_front(oldEdge2);
                }
                if (find(newEdges.begin(), newEdges.end(), newEdge) == newEdges.end()) {
                        newEdges.push_front(newEdge);
                }
        }

        if (IsDeletableNode(node) && IsConsistent(oldEdges, newEdges)){
                DeleteTwoDegreeNodeInconsistent(node);
                FixNodeIDs();

                assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");

                return true;
        }
        else {
                return false;
        }
}

bool EditableMesh::DeleteTwoDegreeNodeInconsistent(Node* node)
{
        //Update owners of node
        list<NeighborNodes> neighbors = m_mesh->GetNodeOwningNeighbors(node);

        if (neighbors.size() != 2) {
                return false; //Tried to delete a node that isn't connected to exactly two other nodes.
        }

        auto nodeOwningWalls = m_mesh->GetNodeOwningWalls(node);
        if (nodeOwningWalls.size() == 1) {
                nodeOwningWalls.front()->SetLengthDirty();
        }
        else {
                assert(nodeOwningWalls.size() == 0 && m_walls.size() == 0
                        && "A two degree node with no wall owners is only possible for one cell in the mesh (no walls).");
        }

        Node* n1 = neighbors.front().GetNb1();
        Node* n2 = neighbors.front().GetNb2();
        NeighborNodes nb1 = neighbors.front();
        NeighborNodes nb2 = neighbors.back();

        nb1.GetCell()->GetNodes().erase(std::find(nb1.GetCell()->GetNodes().begin(), nb1.GetCell()->GetNodes().end(), node));
        nb1.GetCell()->SetGeoDirty();
        nb2.GetCell()->GetNodes().erase(std::find(nb2.GetCell()->GetNodes().begin(), nb2.GetCell()->GetNodes().end(), node));
        nb2.GetCell()->SetGeoDirty();

        //Lambda to test whether the given neighbor contains 'node'.
        auto is_nb = [node](const NeighborNodes& nb) {
                return (nb.GetNb1() == node || nb.GetNb2() == node);
        };

        //Update owners of the first neighbor of node.
        {
                list<NeighborNodes>& owners = m_mesh->GetNodeOwningNeighbors(n1);
                auto owner = find_if(owners.begin(), owners.end(), is_nb);
                while (owner != owners.end()) {
                        if (owner->GetNb1() == node) {
                                const NeighborNodes tmp(owner->GetCell(), n2, owner->GetNb2());
                                *owner = tmp;
                        }
                        else if (owner->GetNb2() == node) {
                                const NeighborNodes tmp(owner->GetCell(), owner->GetNb1(), n2);
                                *owner = tmp;
                        }
                        owner = find_if(owner, owners.end(), is_nb);
                }
        }

        //Update owners of the second neighbor of node.
        {
                list<NeighborNodes>& owners = m_mesh->GetNodeOwningNeighbors(n2);
                auto owner = find_if(owners.begin(), owners.end(), is_nb);
                while (owner != owners.end()) {
                        if (owner->GetNb1() == node) {
                                const NeighborNodes tmp(owner->GetCell(), n1, owner->GetNb2());
                                *owner = tmp;
                        }
                        else if (owner->GetNb2() == node) {
                                const NeighborNodes tmp(owner->GetCell(), owner->GetNb1(), n1);
                                *owner = tmp;
                        }
                        owner = find_if(owner, owners.end(), is_nb);
                }
        }

        m_mesh->RemoveNodeFromOwnerLists(node);
        auto it = m_nodes.begin();
        while (it != m_nodes.end()) {
                if (node == *it) {
                        break;
                }
                it++;
        }
        m_nodes.erase(it);

        return true;
}

bool EditableMesh::DisplaceNode(Node* node, double x, double y)
{
        double old_x = (*node)[0];
        double old_y = (*node)[1];
        (*node)[0] = x;
        (*node)[1] = y;

        if (IsConsistent(node)) {
                for (const auto &nb : m_mesh->GetNodeOwningNeighbors(node)) {
                        nb.GetCell()->SetGeoDirty();
                }
                for (const auto &wall : m_mesh->GetNodeOwningWalls(node)) {
                        wall->SetLengthDirty();
                }

                return true;
        }
        else {
                (*node)[0] = old_x;
                (*node)[1] = old_y;
                return false;
        }
}

void EditableMesh::FixCellIDs()
{
        auto cmp = [](Cell* c1, Cell* c2) {
                return c1->GetIndex() < c2->GetIndex();
        };

        sort(m_cells.begin(), m_cells.end(), cmp);
        for (unsigned int i = 0; i < m_cells.size(); i++) {
                m_cells[i]->m_index = i;
        }

        m_mesh->m_cell_next_id = m_cells.size();
}

void EditableMesh::FixNodeIDs()
{
        auto cmp = [](Node* n1, Node* n2) {
                return n1->GetIndex() < n2->GetIndex();
        };

        sort(m_nodes.begin(), m_nodes.end(), cmp);
        for (unsigned int i = 0; i < m_nodes.size(); i++) {
                m_nodes[i]->m_index = i;
        }

        m_mesh->m_node_next_id = m_nodes.size();
}

void EditableMesh::FixWallIDs()
{
        auto cmp = [](Wall* w1, Wall* w2) {
                return w1->GetIndex() < w2->GetIndex();
        };

        sort(m_walls.begin(), m_walls.end(),cmp);
        for (unsigned int i = 0; i < m_walls.size(); i++) {
                m_walls[i]->m_wall_index = i;
        }

        m_mesh->m_wall_next_id = m_walls.size();
}

std::shared_ptr<Mesh>& EditableMesh::GetMesh()
{
        return m_mesh;
}

bool EditableMesh::IsConsistent(const list<QLineF>& excluded, const list<QLineF>& included) const
{
        // Lambda to check whether an edge is in the list of excluded edges.
        auto is_excluded = [excluded](const QLineF& edge) {
                return (find(excluded.begin(), excluded.end(), edge) != excluded.end());
        };

        // Lambda to check whether an edge is in the list of included edges.
        auto is_included = [included](const QLineF& edge) {
                return (find(included.begin(), included.end(), edge) != included.end());
        };

        // Include all the edges in the mesh.
        list<QLineF> edges;
        for (Cell* cell : m_cells) {
                auto cit = make_const_circular(cell->GetNodes());
                do {
                        QLineF newEdge((**cit)[0], (**cit)[1], (**next(cit))[0], (**next(cit))[1]);
                        if (!is_excluded(newEdge) && !is_included(newEdge)
                                && find(edges.begin(), edges.end(), newEdge) == edges.end()) {
                                edges.push_back(newEdge);
                        }
                } while (++cit != cell->GetNodes().begin());
        }

        // Check whether any of the included edges intersects with one
        // of the mesh edges, excluding the edges in 'excluded'.
        for (QLineF e1 : edges) {
                for (QLineF e2 : included) {
                        QLineF::IntersectType intersection = e1.intersect(e2, nullptr);
                        if (intersection == QLineF::BoundedIntersection
                                && e1.p1() != e2.p1() && e1.p1() != e2.p2()
                                && e1.p2() != e2.p1() && e1.p2() != e2.p2()) {

                                return false;
                        }
                }
        }

        return true;
}

bool EditableMesh::IsConsistent(Node* node) const
{
        list<Cell*> neighbors;
        for (const NeighborNodes& nb : m_mesh->GetNodeOwningNeighbors(node)) {
                if (nb.GetCell()->GetIndex() != -1 && find(neighbors.begin(), neighbors.end(), nb.GetCell()) == neighbors.end()) {
                        neighbors.push_back(nb.GetCell());
                }
        }

        if (node->IsAtBoundary()) {
                for (Node* n : m_mesh->GetBoundaryPolygon()->GetNodes()) {
                        for (const NeighborNodes& nb : m_mesh->GetNodeOwningNeighbors(n)) {
                                if (nb.GetCell()->GetIndex() != -1 && find(neighbors.begin(), neighbors.end(), nb.GetCell()) == neighbors.end()) {
                                        neighbors.push_back(nb.GetCell());
                                }
                        }
                }
        }

        list<QPolygonF> polygons;
        for (Cell* c : neighbors) {
                QPolygonF polygon;
                for (Node* n : c->GetNodes()) {
                        polygon.append(QPointF((*n)[0], (*n)[1]));
                }
                polygons.push_back(polygon);
        }

        for (auto p1 = polygons.begin(); p1 != polygons.end(); p1++) {
                if (!PolygonUtils::IsSimplePolygon(*p1)) {
                        return false;
                }

                for (auto p2 = next(p1); p2 != polygons.end(); p2++) {
                        if (p1->intersected(*p2).size() >= 3) {
                                return false;
                        }
                }
        }

        return true;
}

bool EditableMesh::IsDeletableCell(Cell* cell) const
{
        if (!cell->HasBoundaryWall() || m_cells.size() <= 1
                || find_if(m_cells.begin(), m_cells.end(), [cell](Cell* c){return c == cell;}) == m_cells.end()) {

                return false;
        }

        for (Wall* wall : cell->GetWalls()) {
                if (wall->IsAtBoundary()) {
                        auto nodes = cell->GetNodes();

                        auto firstNode = (wall->GetC1() == cell ?  wall->GetN2() : wall->GetN1());
                        rotate(nodes.begin(), next(find(nodes.begin(), nodes.end(), firstNode)), nodes.end());
                        auto secondNode = (wall->GetC1() == cell ?  wall->GetN1() : wall->GetN2());

                        for (auto it = nodes.begin(); it != find(nodes.begin(), nodes.end(), secondNode); it++) {
                                if ((*it)->IsAtBoundary()) {
                                        //There's a node not on the wall that's at the boundary.
                                        return false;
                                }
                        }

                        return true;
                }
        }
        assert(false && "There's no boundary wall, while this already has been verified.");

        return false;
}

bool EditableMesh::IsDeletableNode(Node* node) const
{
        const list<NeighborNodes> neighbors = m_mesh->GetNodeOwningNeighbors(node);
        if (neighbors.size() != 2) {
                return false; //Tried to delete a node that isn't connected to exactly two other nodes.
        }

        NeighborNodes nb1 = neighbors.front();
        NeighborNodes nb2 = neighbors.back();

        assert(((nb1.GetNb1() == nb2.GetNb1() && nb1.GetNb2() == nb2.GetNb2())
                || (nb1.GetNb1() == nb2.GetNb2() && nb1.GetNb2() == nb2.GetNb1()))
                && nb1.GetCell() != nb2.GetCell()
                && "Inconsistency in the neighbor list of the given node.");

        if (nb1.GetCell()->GetNodes().size() <= 3 || nb2.GetCell()->GetNodes().size() <= 3) {
                return false; //After deleting the node, it wouldn't be a valid cell anymore.
        }

        return true;
}

void EditableMesh::ReplaceCell(Cell* cell, list<QPolygonF> newCells)
{
        assert(find_if(newCells.begin(), newCells.end(),
                        [](const QPolygonF& p){return p.isClosed();}) == newCells.end()
                && "All the given polygons should be open.");
        assert(find_if(newCells.begin(), newCells.end(),
                        [](const QPolygonF& p){return p.size() < 3;}) == newCells.end()
                && "There's a given polygon which is invalid. (less than three corners)");
        assert(!newCells.empty() && "There are no new cells given.");
        assert(all_of(newCells.begin(), newCells.end(),
                        [](const QPolygonF& p){return !PolygonUtils::IsClockwise(p);})
                && "The given polygons weren't counter-clockwise.");
        assert(all_of(cell->GetNodes().begin(), cell->GetNodes().end(), [newCells](Node* n){
                        for (const QPolygonF& p : newCells) {
                                if (p.contains(QPointF((*n)[0], (*n)[1]))) {
                                        return true;
                                }
                        }
                        return false;
                }) && "Not all the points of the cell were present in the new cells.");

        if (newCells.size() == 1) {
                // newCells contains the old cell.
                return;
        }

        // Lambda to compare two points.
        auto cmp_points = [](const QPointF& p1, const QPointF& p2) {
                return p1 != p2 && (p1.x() < p2.x() || (p1.x() == p2.x() && p1.y() < p2.y()));
        };

        // Lambda to compare two edges.
        auto cmp_edges = [](const Edge& e1, const Edge& e2) {
                return e1.GetFirst()->GetIndex() < e2.GetFirst()->GetIndex()
                        || (e1.GetFirst()->GetIndex() == e2.GetFirst()->GetIndex()
                                && e1.GetSecond()->GetIndex() < e2.GetSecond()->GetIndex());
        };

        // Lambda to check whether to points are equal with a small perturbation.
        auto is_perturbed = [](const QPointF& p1, const QPointF& p2) {
                return fabs(p1.x() - p2.x()) < g_accuracy && fabs(p1.y() - p2.y()) < g_accuracy;
        };

        // Lambda to perturb a line.
        auto perturb_line = [](QLineF line) {
                line.setLength(line.length() + g_accuracy);
                line.setPoints(line.p2(), line.p1());
                line.setLength(line.length() + g_accuracy);
                line.setPoints(line.p2(), line.p1());
                return line;
        };

        // Delete the old cell and keep track of its nodes and its edge owners.
        list<Cell*> cells;
        map<QPointF, Node*, decltype(cmp_points)> newNodes(cmp_points);
        map<Edge, Cell*, decltype(cmp_edges)> edgeOwners(cmp_edges);


        auto edge_cit = make_const_circular(cell->GetNodes());
        do {
                Edge edge(*edge_cit, *next(edge_cit));
                edgeOwners[edge] = m_mesh->FindEdgeNeighbor(cell, edge);
        } while (++edge_cit != cell->GetNodes().begin());

        list<QPointF> points;
        for (Node* node : cell->GetNodes()) {
                points.push_back(QPointF((*node)[0], (*node)[1]));
                newNodes[QPointF((*node)[0], (*node)[1])] = node;
        }

        for (const QPolygonF& polygon : newCells) {
                vector<unsigned int> nodeIDs;
                for (const QPointF& point : polygon) {
                        points.push_back(point);

                        if (newNodes.find(point) == newNodes.end()) {
                                newNodes[point] = m_mesh->BuildNode({{point.x(), point.y(), 0}}, false);
                        }
                        nodeIDs.push_back(newNodes[point]->GetIndex());
                }

                Cell* c = m_mesh->BuildCell(nodeIDs);
                c->CellAttributes::operator=(*cell);
                cells.push_back(c);
        }

        std::list<Wall*> newWalls;
        Node* endpoint = nullptr;       //Only used for the special case where we have no walls a priori.
        for (const QPolygonF& polygon : newCells) {
                auto polygonStd = polygon.toStdVector();
                auto polygon_cit = make_const_circular(polygonStd);
                do {
                        auto cell_cit = make_circular(cell->GetNodes());
                        do {
                                QPointF cell_p1((**cell_cit)[0], (**cell_cit)[1]);
                                QPointF cell_p2((**next(cell_cit))[0], (**next(cell_cit))[1]);
                                QLineF cell_line = perturb_line(QLineF(cell_p1, cell_p2));

                                QLineF polygon_line = perturb_line(QLineF(*polygon_cit, *next(polygon_cit)));

                                QPointF intersection;
                                if (polygon_line.intersect(cell_line, &intersection) == QLineF::BoundedIntersection) {
                                        if (is_perturbed(intersection, *polygon_cit)) {
                                                intersection = *polygon_cit;
                                                assert(newNodes.find(intersection) != newNodes.end()
                                                        && "All intersections should already have been added.");
                                        }
                                        else if (is_perturbed(intersection, *next(polygon_cit))) {
                                                intersection = *next(polygon_cit);
                                                assert(newNodes.find(intersection) != newNodes.end()
                                                        && "All intersections should already have been added.");
                                        }
                                        else { //intersection with perturbed, parallel line in the middle.
                                                continue;
                                        }

                                        if (intersection != cell_p1 && intersection != cell_p2) {
                                                Edge edge(*cell_cit, *next(cell_cit));
                                                Edge edge1(*cell_cit, newNodes[intersection]);
                                                Edge edge2(newNodes[intersection], *next(cell_cit));
                                                Cell* neighbor = edgeOwners[edge];

                                                if (m_mesh->GetWalls().size() == 0 && endpoint == nullptr) {
                                                        endpoint = newNodes[intersection];
                                                }
                                                else if (m_mesh->GetWalls().size() == 0) {
                                                        Wall* newWall = m_mesh->BuildWall(newNodes[intersection], endpoint, cell, m_mesh->GetBoundaryPolygon());
                                                        cell->AddWall(newWall);
                                                        neighbor->AddWall(newWall);
                                                        newWalls.push_back(newWall);
                                                        newWall = m_mesh->BuildWall(endpoint, newNodes[intersection], cell, m_mesh->GetBoundaryPolygon());
                                                        cell->AddWall(newWall);
                                                        neighbor->AddWall(newWall);
                                                        newWalls.push_back(newWall);
                                                }
                                                else {
                                                        auto wall = find_if(cell->GetWalls().begin(), cell->GetWalls().end(), [edge](Wall* w){return w->IncludesEdge(edge);});
                                                        assert(wall != cell->GetWalls().end()
                                                                && find_if(next(wall), cell->GetWalls().end(), [edge](Wall* w){return w->IncludesEdge(edge);}) == cell->GetWalls().end()
                                                                && "There should be exactly one wall associated with this edge.");

                                                        // We know in which order the endpoints should be, as the edges are defined in the same order as the nodes in the deleted cell.
                                                        Wall* newWall;
                                                        if ((*wall)->GetC1() == cell) {
                                                                newWall = m_mesh->BuildWall(newNodes[intersection], (*wall)->GetN2(), cell, neighbor);
                                                        }
                                                        else { // (*wall)->GetC2() == cell
                                                                newWall = m_mesh->BuildWall(newNodes[intersection], (*wall)->GetN2(), neighbor, cell);
                                                        }
                                                        (*wall)->SetN2(newNodes[intersection]);
                                                        (*wall)->SetLengthDirty();
                                                        cell->AddWall(newWall);
                                                        neighbor->AddWall(newWall);

                                                        newWalls.push_back(*wall);
                                                        newWalls.push_back(newWall);
                                                }

                                                edgeOwners[edge1] = neighbor;
                                                edgeOwners[edge2] = neighbor;
                                                edgeOwners.erase(edge);

                                                if (neighbor->GetIndex() != -1) {
                                                        m_mesh->AddNodeToCell(neighbor, newNodes[intersection], *next(cell_cit), *cell_cit);
                                                }
                                                else {
                                                        m_mesh->AddNodeToCell(neighbor, newNodes[intersection], *cell_cit, *next(cell_cit));
                                                }
                                                m_mesh->AddNodeToCell(cell, newNodes[intersection], *cell_cit, *next(cell_cit));

                                                cell_cit = cell->GetNodes().begin();
                                        }
                                        else if (newNodes.find(*polygon_cit) != newNodes.end() && newNodes.find(*next(polygon_cit)) != newNodes.end()) {
                                                Edge polygon_edge(newNodes[*polygon_cit], newNodes[*next(polygon_cit)]);

                                                if (!cell->HasEdge(polygon_edge)) {
                                                        //When an edge of the new polygon goes right through a node that already existed in the old cell (and the edge isn't part of the old cell).

                                                        Edge edge(*cell_cit, *next(cell_cit));
                                                        Cell* neighbor = edgeOwners[edge];

                                                        if (m_mesh->GetWalls().size() == 0 && endpoint == nullptr) {
                                                                endpoint = newNodes[intersection];
                                                        }
                                                        else if (m_mesh->GetWalls().size() == 0) {
                                                                if (endpoint != newNodes[intersection]) {
                                                                        Wall* newWall = m_mesh->BuildWall(newNodes[intersection], endpoint, cell, m_mesh->GetBoundaryPolygon());
                                                                        cell->AddWall(newWall);
                                                                        neighbor->AddWall(newWall);
                                                                        newWalls.push_back(newWall);
                                                                        newWall = m_mesh->BuildWall(endpoint, newNodes[intersection], cell, m_mesh->GetBoundaryPolygon());
                                                                        cell->AddWall(newWall);
                                                                        neighbor->AddWall(newWall);
                                                                        newWalls.push_back(newWall);
                                                                }
                                                        }
                                                        else {
                                                                auto wall = find_if(cell->GetWalls().begin(), cell->GetWalls().end(), [edge](Wall* w){return w->IncludesEdge(edge);});
                                                                assert(wall != cell->GetWalls().end()
                                                                        && find_if(next(wall), cell->GetWalls().end(), [edge](Wall* w){return w->IncludesEdge(edge);}) == cell->GetWalls().end()
                                                                        && "There should be exactly one wall associated with this edge.");

                                                                if ((*wall)->GetN1() != newNodes[intersection] && (*wall)->GetN2() != newNodes[intersection]) {
                                                                        Wall* newWall;
                                                                        if ((*wall)->GetC1() == cell) {
                                                                                newWall = m_mesh->BuildWall(newNodes[intersection], (*wall)->GetN2(), cell, neighbor);
                                                                        }
                                                                        else { // (*wall)->GetC2() == cell
                                                                                newWall = m_mesh->BuildWall(newNodes[intersection], (*wall)->GetN2(), neighbor, cell);
                                                                        }
                                                                        (*wall)->SetN2(newNodes[intersection]);
                                                                        (*wall)->SetLengthDirty();
                                                                        cell->AddWall(newWall);
                                                                        neighbor->AddWall(newWall);

                                                                        newWalls.push_back(*wall);
                                                                        newWalls.push_back(newWall);
                                                                }
                                                        }
                                                }
                                        }
                                }
                        } while (++cell_cit != cell->GetNodes().begin());
                } while (++polygon_cit != polygonStd.begin());
        }

        // Update neighbors.
        for (const auto& n : newNodes) {
                m_mesh->GetNodeOwningNeighbors(n.second).remove_if(
                        [cell](const NeighborNodes& nb){return nb.GetCell() == cell;});
        }

        // Update wall owners.
        for (const auto& w : newWalls) {
                m_mesh->UpdateNodeOwningWalls(w);
        }

        // Create new walls in between new cells or reassign the walls from the old cell to the correct cells.
        map<Edge, Cell*, decltype(cmp_edges)> neighbors(cmp_edges);
        for (Cell* c : cells) {
                auto c_cit = make_const_circular(c->GetNodes());
                do {
                        Edge edge(*c_cit, *next(c_cit));
                        Edge reversedEdge(*next(c_cit), *c_cit);

                        auto wall = find_if(cell->GetWalls().begin(), cell->GetWalls().end(),
                                                [edge](Wall* w){return w->IncludesEdge(edge);});
                        if (wall != cell->GetWalls().end()) {
                                assert(find_if(next(wall), cell->GetWalls().end(),
                                                [edge](Wall* w){return w->IncludesEdge(edge);})
                                                == cell->GetWalls().end()
                                        && "There should be exactly one wall associated with this edge.");

                                if ((*wall)->GetC1() == cell) {
                                        (*wall)->SetC1(c);
                                }
                                else {
                                        (*wall)->SetC2(c);
                                }
                                cell->ReassignWall(*wall, c);
                                m_mesh->UpdateNodeOwningWalls(*wall);
                        }
                        else if (edgeOwners.find(edge) == edgeOwners.end()
                                        && edgeOwners.find(reversedEdge) == edgeOwners.end()) {
                                if (neighbors.find(reversedEdge) == neighbors.end()
                                                && neighbors.find(edge) == neighbors.end()) {
                                        neighbors[edge] = c;
                                }
                                else {
                                        Edge owner = (neighbors.find(edge) != neighbors.end()
                                                ? edge : reversedEdge);

                                        Cell* neighbor = neighbors[owner];
                                        neighbors.erase(owner);

                                        Wall* newWall = m_mesh->BuildWall(owner.GetFirst(),
                                                owner.GetSecond(), neighbor, c);

                                        c->AddWall(newWall);
                                        neighbor->AddWall(newWall);
                                        m_mesh->UpdateNodeOwningWalls(newWall);
                                }
                        }
                } while (++c_cit != c->GetNodes().begin());
        }

        m_mesh->RemoveFromNeighborList(cell);

        // Set the new nodes at the boundary if they should be.
        for (const auto& n : newNodes) {
                SetAtBoundary(n.second);
        }

        // Construct the neighbor lists.
        for (Cell* c : cells) {
                m_mesh->ConstructNeighborList(c);
        }

        list<Cell*> neighborCells;
        for (const auto& owner : edgeOwners) {
                auto neighbor = find(neighborCells.begin(), neighborCells.end(), owner.second);
                if (neighbor == neighborCells.end()) {
                        m_mesh->ConstructNeighborList(owner.second);
                        neighborCells.push_front(owner.second);
                }
        }

        auto meshCell = find_if(m_mesh->GetCells().begin(), m_mesh->GetCells().end(),
                                [cell](Cell* c){return c == cell;});
        m_mesh->GetCells().erase(meshCell);
        FixCellIDs();

        assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");
}

void EditableMesh::SetAtBoundary(Node* node)
{
        auto boundary = find_if(
                m_mesh->GetNodeOwningNeighbors(node).begin(),
                m_mesh->GetNodeOwningNeighbors(node).end(),
                [](const NeighborNodes& nb){return nb.GetCell()->GetIndex() == -1;});
        node->SetAtBoundary(boundary != m_mesh->GetNodeOwningNeighbors(node).end());
}

Node* EditableMesh::SplitEdge(const Edge& edge)
{
        assert(m_mesh->GetEdgeOwners(edge).size() > 0 && "Edge doesn't exist in the mesh.");

        Wall* wall = m_mesh->FindWall(edge);
        Node* node = m_mesh->BuildNode({{((*edge.GetFirst())[0] + (*edge.GetSecond())[0]) / 2, ((*edge.GetFirst())[1] + (*edge.GetSecond())[1]) / 2, 0}}, m_mesh->IsAtBoundary(&edge));
        for (NeighborNodes neighbor : m_mesh->GetEdgeOwners(edge)) {
                if (find(neighbor.GetCell()->GetNodes().begin(), neighbor.GetCell()->GetNodes().end(), node) == neighbor.GetCell()->GetNodes().end()) {
                        m_mesh->AddNodeToCell(neighbor.GetCell(), node, edge.GetFirst(), edge.GetSecond());
                }
        }

        if (m_mesh->GetWalls().size() > 0) {
                m_mesh->UpdateNodeOwningWalls(wall);
        }

        assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");
        return node;
}

vector<Cell*> EditableMesh::SplitCell(Node* node1, Node* node2)
{
        QLineF cut((*node1)[0], (*node1)[1], (*node2)[0], (*node2)[1]);

        list<Cell*> owners;
        for (auto nb : m_mesh->GetNodeOwningNeighbors(node1)) {
                if (nb.GetCell()->GetIndex() != -1) {
                        owners.push_back(nb.GetCell());
                }
        }
        for (auto nb : m_mesh->GetNodeOwningNeighbors(node2)) {
                if (nb.GetCell()->GetIndex() != -1) {
                        owners.push_back(nb.GetCell());
                }
        }
        owners.sort([](Cell* c1, Cell* c2) {return c1 < c2;});

        auto eql = [](Cell* c1, Cell* c2) {
                return c1 == c2;
        };

        list<Cell*> commonOwners;
        auto owner = adjacent_find(owners.begin(), owners.end(), eql);
        while (owner != owners.end()) {
                commonOwners.push_back(*owner);
                owner = adjacent_find(++owner, owners.end(), eql);
        }

        if (commonOwners.empty()) {
                return vector<Cell*>();
        }

        list<Cell*> possibilities;
        for (Cell* c : commonOwners) {
                QPolygonF polygon;
                for (Node* n : c->GetNodes()) {
                        polygon.push_back(QPointF((*n)[0], (*n)[1]));
                }
                //if (polygon.containsPoint(center, Qt::OddEvenFill)) {
                auto polygons = PolygonUtils::SlicePolygon(polygon, cut);
                if (polygons.size() == 2) {
                        possibilities.push_back(c);
                }
        }
        assert(possibilities.size() == 1 && "There can only be one cell that can be split.");

        auto cell = find_if(m_mesh->GetCells().begin(), m_mesh->GetCells().end(),
                [possibilities](Cell* c){return c == possibilities.front();});

        return SplitCell(*cell, node1, node2);
}

vector<Cell*> EditableMesh::SplitCell(Cell* cell, Node* node1, Node* node2)
{
        assert(cell != nullptr && node1 != nullptr
                && node2 != nullptr && "The input isn't properly initialized.");

        vector<Cell*> newCells;

        // Input verification.
        if (!CanSplitCell(node1, node2, cell)) {
                return newCells;
        }
        else {
                for (NeighborNodes nb : m_mesh->GetNodeOwningNeighbors(node1)) {
                        if (nb.GetNb1() == node2 || nb.GetNb2() == node2) {
                                newCells.push_back(cell);
                                return newCells;
                        }
                }

                // Check whether the line between node1 and node2 intersects with
                // one of the edges of the cell (except for the endpoints).
                // If this is the case, return the old cell.
                QLineF cut((*node1)[0], (*node1)[1], (*node2)[0], (*node2)[1]);
                auto cit = make_const_circular(cell->GetNodes());
                do {
                        QPointF point1((**cit)[0], (**cit)[1]);
                        QPointF point2((**next(cit))[0], (**next(cit))[1]);
                        QPointF intersection;
                        if (cut.intersect(QLineF(point1, point2), &intersection) == QLineF::BoundedIntersection
                                && intersection != point1 && intersection != point2) {
                                newCells.push_back(cell);
                                return newCells;
                        }
                } while (++cit != cell->GetNodes().begin());
        }

        CellDivider cellDivider(m_mesh.get());
        for (Cell* c : cellDivider.GeometricSplitCell(cell, node1, node2)) {
                auto newCell = find_if(m_mesh->GetCells().begin(), m_mesh->GetCells().end(),
                                                [c](Cell* c2){ return c2 == c; });
                assert(newCell != m_mesh->GetCells().end() && "The new cell can't be found in the mesh.");
                newCells.push_back(*newCell);
        }

        assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");
        return newCells;
}

} // namespace