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

#include "Cell.h"
#include "BoundaryType.h"
#include "Edge.h"
#include "MeshState.h"
#include "Node.h"
#include "WallType.h"

#include "util/container/circular_iterator.h"
#include "util/misc/Exception.h"
#include "util/misc/log_debug.h"

#include <boost/property_tree/ptree.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/optional.hpp>

#include <algorithm>
#include <array>
#include <cmath>
#include <functional>
#include <iomanip>
#include <ios>
#include <iostream>
#include <iterator>
#include <list>
#include <memory>
#include <set>
#include <sstream>
#include <string>
#include <tuple>
#include <vector>

using namespace boost::property_tree;
using boost::optional;
using boost::lexical_cast;
using namespace std;
using namespace SimPT_Sim;
using namespace SimPT_Sim::Container;
using namespace SimPT_Sim::Util;

namespace SimPT_Sim {

namespace {
        struct null_deleter {
                void operator()(const void*) {}
        };
}

Mesh::Mesh(unsigned int num_chemicals)
        :  m_boundary_polygon(nullptr),
           m_cell_next_id(0),
           m_node_next_id(0),
           m_num_chemicals(num_chemicals),
           m_wall_next_id(0)
{
}

void Mesh::AddConsecutiveNodeToCell(Cell* c, Node* n, Node* nb1, Node* nb2)
{
        c->GetNodes().push_back(n);
        GetNodeOwningNeighbors(n).push_back(NeighborNodes(c, nb1, nb2));
}

void Mesh::AddNodeToCell(Cell* c, Node* n, Node* nb1, Node* nb2)
{
        assert(std::find(c->GetNodes().begin(), c->GetNodes().end(), nb1) != c->GetNodes().end() && "Neighbor 1 not present in cell!");
        assert(std::find(c->GetNodes().begin(), c->GetNodes().end(), nb2) != c->GetNodes().end() && "Neighbor 2 not present in cell!");

        auto ins_pos = find(c->GetNodes().begin(), c->GetNodes().end(), nb1);

        // ... second node comes before or after it ...
        if (*(++ins_pos != c->GetNodes().end() ? ins_pos : c->GetNodes().begin()) != nb2) {
                c->GetNodes().insert(
                                ((ins_pos--) != c->GetNodes().begin() ?
                                                ins_pos : (--c->GetNodes().end())),
                                n);
                //cells[c->cell].nodes.insert(--ins_pos, new_node->index);
                // .. set the neighbors of the new node ...
                // in this case e.second and e.first are inverted
                GetNodeOwningNeighbors(n).push_back(
                                NeighborNodes(c, nb2, nb1));
        } else {
                // insert before second node, so leave ins_pos as it is, that is incremented
                c->GetNodes().insert(ins_pos, n);
                // .. set the neighbors of the new node ...
                GetNodeOwningNeighbors(n).push_back(
                                NeighborNodes(c, nb1, nb2));
        }

        // Lambda to test cell index equality
        auto is_nb = [c](const NeighborNodes& n) {
                return (n.GetCell()->GetIndex() == c->GetIndex());
        };

        // Find cell c in NeighborNodes owners of Node edge.first and correct the record
        {
                auto& own = GetNodeOwningNeighbors(nb1);
                const auto cpos = find_if(own.begin(), own.end(), is_nb);
                if (cpos->GetNb1() == nb2) {
                        const NeighborNodes tmp(cpos->GetCell(), n, cpos->GetNb2());
                        *cpos = tmp;
                } else if (cpos->GetNb2() == nb2) {
                        const NeighborNodes tmp(cpos->GetCell(), cpos->GetNb1(), n);
                        *cpos = tmp;
                }
        }

        // Same for Node edge.second
        {
                auto& own = GetNodeOwningNeighbors(nb2);
                const auto cpos = find_if(own.begin(), own.end(), is_nb);
                if (cpos->GetNb1() == nb1) {
                        const NeighborNodes tmp(cpos->GetCell(), n, cpos->GetNb2());
                        *cpos = tmp;
                } else if (cpos->GetNb2() == nb1) {
                        const NeighborNodes tmp(cpos->GetCell(), cpos->GetNb1(), n);
                        *cpos = tmp;
                }
        }
}

Cell* Mesh::BuildBoundaryPolygon(const vector<unsigned int>& node_ids)
{
        auto new_cell = m_cell_storage.emplace_back(-1, GetNumChemicals());
        //m_boundary_polygon = Cell*(new_cell, null_deleter());
        m_boundary_polygon = new_cell;

        unsigned int const nnodes = node_ids.size();
        for (unsigned int i = 0; i < nnodes; i++) {
                assert((GetNodes()[node_ids[i]] != nullptr) && "Node has to exist already!");
                assert((GetNodes()[node_ids[(nnodes + i - 1) % nnodes]] != nullptr) && "Node has to exist already!");
                assert((GetNodes()[node_ids[(i + 1) % nnodes]] != nullptr) && "Node has to exist already!");
                AddConsecutiveNodeToCell(m_boundary_polygon, GetNodes()[node_ids[i]],
                                GetNodes()[node_ids[(nnodes + i - 1) % nnodes]],
                                GetNodes()[node_ids[(i + 1) % nnodes]]);
        }

        return m_boundary_polygon;
}

Cell* Mesh::BuildCell()
{
        auto new_cell = m_cell_storage.emplace_back(m_cell_next_id++, GetNumChemicals());
        //auto c = Cell*(new_cell, null_deleter());
        auto c = new_cell;
        m_cells.push_back(c);
        return c;
}

Cell* Mesh::BuildCell(const std::vector<unsigned int>& node_ids)
{
        return BuildCell(m_cell_next_id, node_ids);
}

Cell* Mesh::BuildCell(unsigned int id, const vector<unsigned int>& node_ids)
{
        assert((id >= m_cell_next_id) && "Cell id not increasing!");
        auto new_cell = m_cell_storage.emplace_back(id++, GetNumChemicals());
        //auto c = Cell*(new_cell, null_deleter());
        auto c = new_cell;
        m_cell_next_id = max(m_cell_next_id, id);
        m_cells.push_back(c);

        unsigned int const nnodes = node_ids.size();
        for (unsigned int i = 0; i < nnodes; i++) {
                assert((GetNodes()[node_ids[i]] != nullptr) && "Node has to exist already!");
                assert((GetNodes()[node_ids[(nnodes + i - 1) % nnodes]] != nullptr) && "Node has to exist already!");
                assert((GetNodes()[node_ids[(i + 1) % nnodes]] != nullptr) && "Node has to exist already!");
                AddConsecutiveNodeToCell(c, GetNodes()[node_ids[i]],
                                GetNodes()[node_ids[(nnodes + i - 1) % nnodes]],
                                GetNodes()[node_ids[(i + 1) % nnodes]]);
        }

        return c;
}

Node* Mesh::BuildNode(const array<double, 3>& vec, bool at_boundary)
{
        auto new_node = m_node_storage.emplace_back(m_node_next_id++, vec, at_boundary);
        //auto n = Node*(new_node, null_deleter());
        auto n = new_node;
        m_nodes.push_back(n);
        return n;
}

Node* Mesh::BuildNode(unsigned int id, const array<double, 3>& vec, bool at_boundary)
{
        assert((id >= m_node_next_id) && "Node id not increasing!");
        auto new_node = m_node_storage.emplace_back(id++, vec, at_boundary);
        //auto n = Node*(new_node, null_deleter());
        auto n = new_node;
        m_node_next_id = max(m_node_next_id, id);
        m_nodes.push_back(n);
        return n;
}

Wall* Mesh::BuildWall(Node* n1, Node* n2, Cell* c1, Cell* c2)
{
        assert( (n1 != nullptr) && "Wall::Wall> Attempting to build a wall with nullptr node!" );
        assert( (n2 != nullptr) && "Wall::Wall> Attempting to build a wall with nullptr node!" );
        assert( (c1 != nullptr) && "Wall::Wall> Attempting to build a wall with nullptr cell!" );
        assert( (c2 != nullptr) && "Wall::Wall> Attempting to build a wall with nullptr cell!" );
        assert( (n1 != n2) && "Wall::Wall> Attempting to build a wall between identical nodes!" );
        assert( (c1 != c2) && "Wall::Wall> Attempting to build a wall between identical cells!" );
        //if ( m_model_name != "Wortel" ){};//DDV
        auto new_wall = m_wall_storage.emplace_back(m_wall_next_id++, n1, n2, c1, c2, GetNumChemicals()/*,isWortel*/);
        //auto wall = Wall*(new_wall, null_deleter());
        auto wall = new_wall;
        m_walls.push_back(wall);

        return wall;
}

Wall* Mesh::BuildWall(unsigned int id, Node* n1, Node* n2, Cell* c1, Cell* c2)
{
        assert( (n1 != nullptr) && "Wall::Wall> Attempting to build a wall with nullptr node!" );
        assert( (n2 != nullptr) && "Wall::Wall> Attempting to build a wall with nullptr node!" );
        assert( (c1 != nullptr) && "Wall::Wall> Attempting to build a wall with nullptr cell!" );
        assert( (c2 != nullptr) && "Wall::Wall> Attempting to build a wall with nullptr cell!" );
        assert( (n1 != n2) && "Wall::Wall> Attempting to build a wall between identical nodes!" );
        assert( (c1 != c2) && "Wall::Wall> Attempting to build a wall between identical cells!" );
        //As above...
        auto new_wall = m_wall_storage.emplace_back(id++, n1, n2, c1, c2, GetNumChemicals()/*,*/);
        //auto wall = Wall*(new_wall, null_deleter());
        auto wall = new_wall;
        m_wall_next_id = max(m_wall_next_id, id);
        m_walls.push_back(wall);

        return wall;
}

void Mesh::ConstructNeighborList(Cell* cell)
{
        list<Cell*>& neighbors = m_wall_neighbors[cell];
        neighbors.clear();

        const list<Wall*>& walls = cell->GetWalls();
        for (auto const& wall : walls) {
                if (wall->GetC1() != cell && !wall->GetC1()->IsBoundaryPolygon()) {
                        neighbors.push_back(wall->GetC1());
                } else if (wall->GetC2() != cell && !wall->GetC2()->IsBoundaryPolygon()) {
                        neighbors.push_back(wall->GetC2());
                }
        }
}

Cell* Mesh::FindEdgeNeighbor(Cell* cell, Edge edge) const
{
        Cell* neighbor_cell = nullptr;

        // Push all cells owning the two nodes of the divides edges onto a list
        list<NeighborNodes> owners;
        auto& owners_first = GetNodeOwningNeighbors(edge.GetFirst());
        copy(owners_first.begin(), owners_first.end(), back_inserter(owners));
        auto& owners_second = GetNodeOwningNeighbors(edge.GetSecond());
        copy(owners_second.begin(), owners_second.end(), back_inserter(owners));

        // Lambdas to compare neighbors.
        auto cmp = [](const NeighborNodes& n1, const NeighborNodes& n2) {
                return n1.GetCell()->GetIndex() < n2.GetCell()->GetIndex();
        };
        auto eql = [](const NeighborNodes& n1, const NeighborNodes& n2) {
                return n1.GetCell()->GetIndex() == n2.GetCell()->GetIndex();
        };

        // Find first non-self duplicate in the owners:
        // cells owning the same two nodes share an edge with me
        owners.sort(cmp);
        list<NeighborNodes>::iterator it;
        for (it = owners.begin(); it != owners.end(); ++it) {
                it = adjacent_find(it, owners.end(), eql);
                if (it == owners.end()) {
                        break; // bail if reach the end of the list
                }
                if (it->GetCell()->GetIndex() != cell->GetIndex()) {
                        // This is necessary because two endpoints can have more than two similar cells neighbors.
                        // For example:
                        //
                        //      1---2---3---4
                        //      |   |___|   |
                        //      |___________|
                        //
                        //      2 and 3 have three cell neighbors the same.

                        bool isCorrectNeighbor = false;

                        auto cit = make_const_circular(it->GetCell()->GetNodes());
                        do {
                                if (Edge(*cit, *next(cit)) == edge) {
                                        isCorrectNeighbor = true;
                                }
                        } while (!isCorrectNeighbor && ++cit != it->GetCell()->GetNodes().begin());

                        if (isCorrectNeighbor) {
                                neighbor_cell = it->GetCell();
                                break;
                        }
                }
        }
        assert((neighbor_cell != nullptr) && "There has to be a neighbor");
        return neighbor_cell;

        //TODO: This doesn't always work (see above why.)

        /*Cell* neighbor_cell = nullptr;

        // Push all cells owning the two nodes of the divides edges onto a list
        list<Neighbor> owners;
        auto& owners_first = GetNodeOwningNeighbors(edge.GetFirst());
        copy(owners_first.begin(), owners_first.end(), back_inserter(owners));
        auto& owners_second = GetNodeOwningNeighbors(edge.GetSecond());
        copy(owners_second.begin(), owners_second.end(), back_inserter(owners));

        // Lambdas to compare neighbors.
        auto cmp = [](const Neighbor& n1, const Neighbor& n2) {
                return n1.GetCell()->GetIndex() < n2.GetCell()->GetIndex();
        };
        auto eql = [](const Neighbor& n1, const Neighbor& n2) {
                return n1.GetCell()->GetIndex() == n2.GetCell()->GetIndex();
        };

        // Find first non-self duplicate in the owners:
        // cells owning the same two nodes share an edge with me
        owners.sort(cmp);
        list<Neighbor>::iterator it;
        for (it = owners.begin(); it != owners.end(); it++) {
                it = adjacent_find(it, owners.end(), eql);
                if (it == owners.end()) {
                        break; // bail if reach the end of the list
                }
                if (it->GetCell()->GetIndex() != cell->GetIndex()) {
                        neighbor_cell = it->GetCell();
                        break;
                }
        }
        assert((neighbor_cell != nullptr) && "There has to be a neighbor");
        return neighbor_cell;
        */
}

Node* Mesh::FindNextBoundaryNode(Node* boundary_node) const
{
        auto current = find(m_boundary_polygon->GetNodes().begin(),
                                m_boundary_polygon->GetNodes().end(), boundary_node);
        if (current != m_boundary_polygon->GetNodes().end()) {
                return (++current != m_boundary_polygon->GetNodes().end())
                        ? *current : *m_boundary_polygon->GetNodes().begin();
        } else {
                return nullptr;
        }
}

Wall* Mesh::FindWall(const Edge& edge) const
{
        if (m_walls.empty()) {
                assert(m_cells.size() == 1 && "The only case without walls in mesh, is single cell.");
                return nullptr;
        }

        auto neighbors = GetEdgeOwners(edge); // neighbors == [C1, C1, C2, C2].
        Cell* cell1 = neighbors.front().GetCell();
        Cell* cell2 = neighbors.back().GetCell();

        //Lambda to check whether the given wall doesn't contain both cell1 and cell2.
        auto not_contain_cells = [cell1, cell2](Wall* w){
                return (w->GetC1() != cell1 && w->GetC2() != cell1)
                        || (w->GetC1() != cell2 && w->GetC2() != cell2);
        };

        auto walls = cell1->GetWalls();
        walls.remove_if(not_contain_cells);
        assert(!walls.empty() && "At least one candidate wall tshould associated with the edge.");

        for (auto w : walls) {
                auto cellNodes = w->GetC1()->GetNodes();
                rotate(cellNodes.begin(), find(cellNodes.begin(),
                        cellNodes.end(), w->GetN1()), cellNodes.end());

                for (auto it = cellNodes.begin(); it != find(cellNodes.begin(), cellNodes.end(), w->GetN2()); it++) {
                        if ((edge.GetFirst() == *it && edge.GetSecond() == *next(it))
                                || (edge.GetSecond() == *it && edge.GetFirst() == *next(it))) {

                                return w;
                        }
                }
        }
        assert(false && "No wall has been found for the given edge.");

        return nullptr;
}

void Mesh::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_cell_next_id = m_cells.size();
}

void Mesh::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_node_next_id = m_nodes.size();
}

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

        sort(m_walls.begin(), m_walls.end(),cmp);
        unsigned int i = 0;
        for (auto wall : m_walls) {
                wall->m_wall_index = i++;
        }

        m_wall_next_id = m_walls.size();
}

vector<NeighborNodes> Mesh::GetEdgeOwners(const Edge& edge) const
{
        list<NeighborNodes> owners;

        // Lambdas  to sort and test-equal NeighborNodes.
        auto cmp = [](const NeighborNodes& n1, const NeighborNodes& n2) {
                return n1.GetCell()->GetIndex() < n2.GetCell()->GetIndex();
        };
        auto eql = [](const NeighborNodes& n1, const NeighborNodes& n2) {
                return n1.GetCell()->GetIndex() == n2.GetCell()->GetIndex();
        };

        // Push all cells owning either of the two nodes of the edge onto a list.
        auto& owners_first = GetNodeOwningNeighbors(edge.GetFirst());
        copy(owners_first.begin(), owners_first.end(), back_inserter(owners));
        auto& owners_second = GetNodeOwningNeighbors(edge.GetSecond());
        copy(owners_second.begin(), owners_second.end(), back_inserter(owners));
        owners.sort(cmp);

        // The duplicates in this list indicate cells owning both nodes.
        vector<NeighborNodes> edge_owners;
        auto n = adjacent_find(owners.begin(), owners.end(), eql);
        while (n != owners.end()) {
                const auto owner_a = *n;
                const auto owner_b = *(++n);
                const Edge e1(owner_a.GetNb1(), owner_b.GetNb2());
                const Edge e2(owner_b.GetNb1(), owner_a.GetNb2());
                if ( e1 == edge || e2 == edge ) {
                        edge_owners.push_back(owner_a);
                        edge_owners.push_back(owner_b);
                }
                n = adjacent_find(n, owners.end(), eql);
        }

        assert(edge_owners.size() == 4 && "An edge can only have two neighboring cells per endpoint.");

        return edge_owners;
}

MeshState Mesh::GetState() const
{
        MeshState state;
        try {
                // Nodes
                state.SetNumNodes(m_nodes.size());
                state.NodeAttributeContainer().Add<int>("fixed");
                state.NodeAttributeContainer().Add<int>("sam");
                state.NodeAttributeContainer().Add<int>("boundary");
                size_t node_idx = 0;
                for (auto const& n : m_nodes) {
                        state.SetNodeID(node_idx, n->GetIndex());
                        state.SetNodeX(node_idx, (*n)[0]);
                        state.SetNodeY(node_idx, (*n)[1]);
                        state.NodeAttributeContainer().Set<int>(node_idx, "fixed", n->IsFixed());
                        state.NodeAttributeContainer().Set<int>(node_idx, "sam", n->IsSam());
                        state.NodeAttributeContainer().Set<int>(node_idx, "boundary", n->IsAtBoundary());
                        ++node_idx;
                }

                // Cells
                state.SetNumCells(m_cells.size());
                state.CellAttributeContainer().Add<int>("type");
                state.CellAttributeContainer().Add<double>("area");
                state.CellAttributeContainer().Add<double>("solute");
                state.CellAttributeContainer().Add<double>("target_area");
                state.CellAttributeContainer().Add<double>("target_length");
                state.CellAttributeContainer().Add<int>("dead");
                state.CellAttributeContainer().Add<int>("div_counter");
                state.CellAttributeContainer().Add<double>("stiffness");
                state.CellAttributeContainer().Add<string>("boundary_type");
                state.CellAttributeContainer().Add<int>("fixed");
                // TODO: a bunch of other attributes: cfr. CellAttributes
                // Cell chemicals are attributes as well
                for (size_t chem_idx = 0; chem_idx < GetNumChemicals(); ++chem_idx) {
                        state.CellAttributeContainer().Add<double>("chem_" + lexical_cast<string>(chem_idx));
                }
                size_t cell_idx = 0;
                for (auto c : m_cells) {
                        // IDs of cells
                        state.SetCellID(cell_idx, c->GetIndex());
                        // Loop over cell's nodes and collect their indices
                        vector<int> cell_nodes;
                        for (auto const& node : c->GetNodes()) {
                                cell_nodes.push_back(node->GetIndex());
                        }
                        state.SetCellNodes(cell_idx, cell_nodes);
                        // Loop over cell's walls and collect their indices
                        vector<int> cell_walls;
                        for (auto const& wall : c->GetWalls()) {
                                cell_walls.push_back(wall->GetIndex());
                        }
                        state.SetCellWalls(cell_idx, cell_walls);
                        // Collect cell attributes
                        state.CellAttributeContainer().Set<int>(
                                cell_idx, "type", c->GetCellType());
                        state.CellAttributeContainer().Set<double>(
                                cell_idx, "area", c->GetArea());
                        state.CellAttributeContainer().Set<double>(
                                cell_idx, "solute", c->GetSolute());
                        state.CellAttributeContainer().Set<double>(
                                cell_idx, "target_area", c->GetTargetArea());
                        state.CellAttributeContainer().Set<double>(
                                cell_idx, "target_length", c->GetTargetLength());
                        state.CellAttributeContainer().Set<int>(
                                cell_idx, "dead", c->IsDead());
                        state.CellAttributeContainer().Set<int>(
                                cell_idx, "div_counter", c->GetDivisionCounter());
                        state.CellAttributeContainer().Set<double>(
                                cell_idx, "stiffness", c->GetStiffness());
                        state.CellAttributeContainer().Set<string>(
                                cell_idx, "boundary_type", ToString(c->GetBoundaryType()));
                        state.CellAttributeContainer().Set<int>(
                                cell_idx, "fixed", c->IsFixed());
                        for (size_t chem_idx = 0; chem_idx < GetNumChemicals(); ++chem_idx) {
                                state.CellAttributeContainer().Set<double>(
                                        cell_idx, "chem_" + lexical_cast<string>(chem_idx),
                                        c->GetChemical(chem_idx));
                        }
                        ++cell_idx;
                }

                // Boundary polygon: is a cell with ID = -1 but is stored separately in MeshState
                vector<int> bp_nodes;
                for (auto const& node : m_boundary_polygon->GetNodes()) {
                        bp_nodes.push_back(node->GetIndex());
                }
                state.SetBoundaryPolygonNodes(bp_nodes);
                vector<int> bp_walls;
                for (auto const& wall : m_boundary_polygon->GetWalls()) {
                        bp_walls.push_back(wall->GetIndex());
                }
                state.SetBoundaryPolygonWalls(bp_walls);

                // The walls
                state.SetNumWalls(m_walls.size());
                state.WallAttributeContainer().Add<double>("length");
                state.WallAttributeContainer().Add<double>("rest_length");
                state.WallAttributeContainer().Add<double>("rest_length_init");
                state.WallAttributeContainer().Add<double>("strength");
                state.WallAttributeContainer().Add<string>("type");
                // Wall transporters are attributes as well
                for (size_t chem_idx = 0; chem_idx < GetNumChemicals(); ++chem_idx) {
                        state.WallAttributeContainer().Add<double>(
                                "trans_1_chem_" + lexical_cast<string>(chem_idx));
                        state.WallAttributeContainer().Add<double>(
                                "trans_2_chem_" + lexical_cast<string>(chem_idx));
                }
                size_t wall_idx = 0;
                for (auto const& w : m_walls) {
                        // IDs of walls
                        state.SetWallID(wall_idx, w->GetIndex());
                        // Connected cells and end nodes
                        state.SetWallNodes(wall_idx, make_pair(w->GetN1()->GetIndex(), w->GetN2()->GetIndex()));
                        state.SetWallCells(wall_idx, make_pair(w->GetC1()->GetIndex(), w->GetC2()->GetIndex()));
                        // Collect wall attributes
                        state.WallAttributeContainer().Set<double>(
                    wall_idx, "length", w->GetLength());
                        state.WallAttributeContainer().Set<double>(
                    wall_idx, "rest_length", w->GetRestLength());
                        state.WallAttributeContainer().Set<double>(
                    wall_idx, "rest_length_init", w->GetRestLengthInit());
                        state.WallAttributeContainer().Set<double>(
                    wall_idx, "strength", w->GetStrength());
                        state.WallAttributeContainer().Set<string>(
                    wall_idx, "type", WallType::ToString(w->GetType()));
                        for (size_t chem_idx = 0; chem_idx < GetNumChemicals(); ++chem_idx) {
                                state.WallAttributeContainer().Set<double>(
                                        wall_idx, "trans_1_chem_" + lexical_cast<string>(chem_idx),
                                        w->GetTransporters1(chem_idx));
                                state.WallAttributeContainer().Set<double>(
                                        wall_idx, "trans_2_chem_" + lexical_cast<string>(chem_idx),
                                        w->GetTransporters2(chem_idx));
                        }
                        ++wall_idx;
                }
        } catch (exception& e) {
                cout << "ERROR in Mesh::GetState(): "  << e.what() << endl;
        }
        return state;
}

bool Mesh::IsAtBoundary(const Edge* e) const
{
        if (e->GetFirst()->IsAtBoundary() && e->GetSecond()->IsAtBoundary()) {
                auto node1 = FindNextBoundaryNode(e->GetFirst());
                auto node2 = FindNextBoundaryNode(e->GetSecond());
                return node1 == e->GetSecond() || node2 == e->GetFirst();
        }
        else {
                return false;
        }
}

bool Mesh::IsAtBoundary(Wall* w) const
{
        bool status = false;
        if (w->GetC1() == m_boundary_polygon || w->GetC2() == m_boundary_polygon) {
                status = true;
        }
        return status;
}

bool Mesh::IsInBoundaryPolygon(Node* n) const
{
        bool status = false;
        const auto& bp_nodes = m_boundary_polygon->GetNodes();
        auto it = find(bp_nodes.begin(), bp_nodes.end(), n);
        if (it != bp_nodes.end()) {
                status = true;
        }
        return status;
}

ostream& Mesh::Print(ostream& os) const {
        static int i = 0;
        i++;
        os << "MESH_PRINT_" << i << endl;
        os << "BOUNDARY_CELL" << endl;
        if (m_boundary_polygon) {
                m_boundary_polygon->Print(os) << endl;
        }
        os << "CELLS" << endl;
        for (auto& cell : m_cells) {
                cell->Print(os) << endl;
        }
        os << "NODES" << endl;
        for (auto& node : m_nodes) {
                node->Print(os) << endl;
        }
        os << "WALLS" << endl;
        for (auto& wall : m_walls) {
                wall->Print(os) << endl;
        }
        os << "NODE_AND_ITS_NEIGHBOR_NODES" << endl;
        {
                set<string> ordered_output;
                ostringstream os;
                for (auto& pair : m_node_owning_neighbors) {
                        os.str("");
                        os << "=====NODE=====" << endl;
                        pair.first->Print(os) << endl;
                        for (auto& n : pair.second) {
                                os << "-----here-comes-an-entry-----" << endl;
                                os << "***CELL: " << endl;
                                n.GetCell()->Print(os) << endl;
                                os << "***NB_1: " << endl;
                                n.GetNb1()->Print(os) << endl;
                                os << "***NB_2: " << endl;
                                os << "-----------------------------" << endl;
                                n.GetNb2()->Print(os) << endl;
                        }
                        os << "==============" << endl;
                        ordered_output.insert(os.str());
                }
                for (string s : ordered_output) {
                        os << s;
                }
        }
        os << "NODE_AND_ITS_WALLS" << endl;
        {
                set<string> ordered_output;
                ostringstream os;
                for (auto& pair : m_node_owning_walls) {
                        os.str("");
                        os << "=====NODE=====" << endl;
                        pair.first->Print(os) << endl;
                        for (auto& wall : pair.second) {
                                os << "-----here-comes-an-entry-----" << endl;
                                os << "***WALL: " << endl;
                                wall->Print(os) << endl;
                                os << "-----------------------------" << endl;
                        }
                        os << "==============" << endl;
                        ordered_output.insert(os.str());
                }
                for (string s : ordered_output) {
                        os << s;
                }
        }
        os << "CELL_AND_ITS_NEIGHBOR_CELLS_HAVING_A_COMMON_WALL" << endl;
        {
                set<string> ordered_output;
                ostringstream os;
                for (auto& pair: m_wall_neighbors) {
                        os.str("");
                        os << "=====CELL=====" << endl;
                        pair.first->Print(os) << endl;
                        for (auto& cell: pair.second) {
                                os << "-----here-comes-an-entry-----" << endl;
                                os << "***NBCELL: " << endl;
                                cell->Print(os) << endl;
                                os << "-----------------------------" << endl;
                        }
                        os << "==============" << endl;
                        ordered_output.insert(os.str());
                }
                for (string s : ordered_output) {
                        os << s;
                }
        }
        return os;
}

void Mesh::RemoveFromNeighborList(Cell* cell)
{
        m_wall_neighbors.erase(cell);
}

void Mesh::RemoveNodeFromOwnerLists(Node* node)
{
        m_node_owning_neighbors.erase(node);
        m_node_owning_walls.erase(node);
}

ptree Mesh::ToPtree() const
{
        // ------------------------- File intro ---------------------------------
        ptree ret_pt;

        // --------------- Cells & BoundaryPolygon ------------------------------
        auto const& cells = GetCells();
        ptree cells_pt;
        cells_pt.put("count", cells.size());
        cells_pt.put("chemical_count", GetNumChemicals());
        for (auto const& cell : cells) {
                cells_pt.add_child("cell_array.cell", cell->ToPtree());
        }
        cells_pt.put_child("boundary_polygon", GetBoundaryPolygon()->ToPtree());
        ret_pt.put_child("cells", cells_pt);

        // --------------------------- Nodes ------------------------------------
        const auto& nodes = GetNodes();
        ptree nodes_pt;
        nodes_pt.put("count", nodes.size());
        for (const auto& node : nodes) {
                nodes_pt.add_child("node_array.node", node->ToPtree());
        }
        ret_pt.put_child("nodes", nodes_pt);

        // --------------------------- Walls ----------------------------------
        const auto& walls = GetWalls();
        ptree walls_pt;
        walls_pt.put("count", walls.size());
        for (auto const& w : walls) {
                walls_pt.add_child("wall_array.wall", w->ToPtree());
        }
        ret_pt.put_child("walls", walls_pt);

        return ret_pt;
}

void Mesh::UpdateNodeOwningNeighbors(Cell* cell)
{
        auto& nodes = cell->GetNodes();
        if (!nodes.empty()) {
                auto it = make_circular(nodes);
                do {
                        auto& owners = GetNodeOwningNeighbors(*it);
                        // Clean existing connections in which this cell participates
                        if (owners.size() > 0) {
                                // remove myself from the neighbor list of the node
                                auto is_nb = [cell](const NeighborNodes& n) {
                                        return (n.GetCell()->GetIndex() == cell->GetIndex());
                                };
                                auto nb_it = find_if(owners.begin(), owners.end(), is_nb);
                                if (nb_it != owners.end()) {
                                        owners.erase(nb_it);
                                }
                        }
                        owners.push_back(NeighborNodes(cell, *prev(it), *next(it)));
                } while (++it != nodes.begin());
        }
}

void Mesh::UpdateNodeOwningWalls(Wall *wall)
{
        auto& nodes = wall->GetC1()->GetNodes();
        auto start = find(nodes.begin(), nodes.end(), wall->GetN1());
        auto stop = find(nodes.begin(), nodes.end(), wall->GetN2());
        assert(start != nodes.end() && "N1 not found in C1 node list");
        assert(stop != nodes.end() && "N2 not found in C1 node list");

        for (auto &endpoint : {start, stop}) {
                auto &walls = m_node_owning_walls[*endpoint];
                walls.remove_if([endpoint] (Wall *owningWall) { return owningWall->GetN1() != *endpoint && owningWall->GetN2() != *endpoint; } );
                if (std::find(walls.begin(), walls.end(), wall) == walls.end()) {
                        walls.push_back(wall);
                }
        }

        auto cit = make_circular(nodes, start);
        while (++cit != stop) {
                m_node_owning_walls[*cit] = { wall }; // Nodes on the wall that are not the endpoints cannot be part of multiple walls
        }
}

}