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

#include "Cell.h"
#include "Mesh.h"
#include "MeshState.h"
#include "Node.h"

#include <boost/algorithm/string/predicate.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/property_tree/ptree.hpp>
#include <boost/property_tree/xml_parser.hpp>
#include <iostream>

namespace SimPT_Sim {

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

shared_ptr<Mesh> MBMBuilder::Build(const MeshState& mesh_state)
{
        unsigned int num_chem = CalculateNumChemicals(mesh_state);
        m_mesh = make_shared<Mesh>(num_chem);

        BuildNodes(mesh_state);
        BuildCells(mesh_state);
        BuildBoundaryPolygon(mesh_state);
        BuildWalls(mesh_state);
        ConnectCellsToWalls(mesh_state);

        m_mesh->ConstructNeighborList(m_mesh->GetBoundaryPolygon());
        m_mesh->UpdateNodeOwningNeighbors(m_mesh->GetBoundaryPolygon());
        for (auto const& cell : m_mesh->GetCells()) {
                m_mesh->ConstructNeighborList(cell);
                m_mesh->UpdateNodeOwningNeighbors(cell);
        }

        // TODO: Nodesets: I don't know what that is or why it would be useful?

        return m_mesh;
}

void MBMBuilder::BuildBoundaryPolygon(const MeshState& mesh_state)
{
        vector<int> bpni = mesh_state.GetBoundaryPolygonNodes();
        vector<unsigned int> bp_nodes_id(bpni.begin(), bpni.end()); // (convert from int IDs to uint IDs)
        auto new_bp = m_mesh->BuildBoundaryPolygon(bp_nodes_id);
        const ptree empty_pt;
        new_bp->Cell::ReadPtree(empty_pt);
}

void MBMBuilder::BuildCells(const MeshState& mesh_state)
{
        unsigned int num_chem = m_mesh->GetNumChemicals();

        // Build cells
        for (size_t cell_idx = 0; cell_idx < mesh_state.GetNumCells(); ++cell_idx) {
                // Construct an empty cell.
                const int cell_id = mesh_state.GetCellID(cell_idx);
                vector<int> cni = mesh_state.GetCellNodes(cell_idx);
                vector<unsigned int> cell_nodes_id(cni.begin(), cni.end()); // (convert from int IDs to uint IDs)
                auto new_cell = m_mesh->BuildCell(cell_id, cell_nodes_id);

                // Cell geometry.
                new_cell->SetGeoDirty();

                // Cell attributes.
                ptree attributes_pt;
                auto cell_attr = mesh_state.CellAttributeContainer();
                attributes_pt.put("boundary",        cell_attr.Get<string>(cell_idx, "boundary_type"));
                attributes_pt.put("dead",            static_cast<bool>(cell_attr.Get<int>(cell_idx, "dead")));
                attributes_pt.put("division_count",  cell_attr.Get<int>(cell_idx, "div_counter"));
                attributes_pt.put("fixed",           cell_attr.Get<int>(cell_idx, "fixed"));
                attributes_pt.put("solute",          cell_attr.Get<double>(cell_idx, "solute"));
                attributes_pt.put("stiffness",       cell_attr.Get<double>(cell_idx, "stiffness"));
                attributes_pt.put("target_area",     cell_attr.Get<double>(cell_idx, "target_area"));
                attributes_pt.put("target_length",   cell_attr.Get<double>(cell_idx, "target_length"));
                attributes_pt.put("type",            cell_attr.Get<int>(cell_idx, "type"));
                for (size_t chem_idx = 0; chem_idx < num_chem; ++chem_idx) {
                        const string chem_name = string("chem_") + boost::lexical_cast<string>(chem_idx);
                        new_cell->SetChemical(chem_idx, cell_attr.Get<double>(cell_idx, chem_name));
                        attributes_pt.add("chemical_array.chemical", cell_attr.Get<double>(cell_idx, chem_name));
                }
                new_cell->CellAttributes::ReadPtree(attributes_pt);
        }
}

void MBMBuilder::BuildNodes(const MeshState& mesh_state)
{
        // Build nodes from nodes contained in mesh_state
        for (size_t node_idx = 0; node_idx < mesh_state.GetNumNodes(); ++node_idx) {
                // Required Node properties & attributes
                const int node_id = mesh_state.GetNodeID(node_idx);
                array<double,3> node_xyz {{mesh_state.GetNodeX(node_idx), mesh_state.GetNodeY(node_idx), 0.0}};
                int node_at_boundary = mesh_state.NodeAttributeContainer().Get<int>(node_idx, "boundary");

                // Basic node building delegated to Mesh::BuildNode
                auto new_node = m_mesh->BuildNode(node_id, node_xyz, node_at_boundary);

                // Required node attributes
                new_node->SetSam(mesh_state.NodeAttributeContainer().Get<int>(node_idx, "sam"));
                new_node->SetFixed(mesh_state.NodeAttributeContainer().Get<int>(node_idx, "fixed"));
        }
}

void MBMBuilder::BuildWalls(const MeshState& mesh_state)
{
        auto& mesh_cells = m_mesh->GetCells();
        auto& mesh_nodes = m_mesh->GetNodes();
        for (size_t wall_idx = 0; wall_idx < mesh_state.GetNumWalls(); ++wall_idx) {
                const int wall_id = mesh_state.GetWallID(wall_idx);
                const pair<int,int> nid = mesh_state.GetWallNodes(wall_idx);
                const pair<int,int> cid = mesh_state.GetWallCells(wall_idx);
                auto n1 = mesh_nodes[nid.first];
                auto n2 = mesh_nodes[nid.second];
                auto c1 = (cid.first  != -1) ? mesh_cells[cid.first]  : m_mesh->GetBoundaryPolygon();
                auto c2 = (cid.second != -1) ? mesh_cells[cid.second] : m_mesh->GetBoundaryPolygon();

                auto new_wall = m_mesh->BuildWall(wall_id, n1, n2, c1, c2);

                // Wall attributes.
                ptree attributes_pt;
                auto wall_attr = mesh_state.WallAttributeContainer();
                attributes_pt.put("rest_length",  wall_attr.Get<double>(wall_idx, "rest_length"));
                attributes_pt.put("rest_length_init",  wall_attr.Get<double>(wall_idx, "rest_length_init"));
                attributes_pt.put("strength",     wall_attr.Get<double>(wall_idx, "strength"));
                attributes_pt.put("type",         wall_attr.Get<string>(wall_idx, "type"));

                // Wall chemicals (transporters)
                // (It is assumed throughout SimPT that walls and cells share the same number of chemicals.
                // Since cells are read first, it is reasonable to rely on correct detection of number of
                // chemicals in ReadCells())
                size_t num_chem = m_mesh->GetNumChemicals();
                for (size_t chem_idx = 0; chem_idx < num_chem; ++chem_idx) {
                        const string trans_1_name = string("trans_1_chem_") + boost::lexical_cast<string>(chem_idx);
                        const string trans_2_name = string("trans_2_chem_") + boost::lexical_cast<string>(chem_idx);
                        new_wall->SetTransporters1(chem_idx, wall_attr.Get<double>(wall_idx, trans_1_name));
                        new_wall->SetTransporters2(chem_idx, wall_attr.Get<double>(wall_idx, trans_2_name));
                        attributes_pt.add("transporter1_array.transporter1", wall_attr.Get<double>(wall_idx, trans_1_name));
                        attributes_pt.add("transporter2_array.transporter2", wall_attr.Get<double>(wall_idx, trans_2_name));
                }
                new_wall->WallAttributes::ReadPtree(attributes_pt);
        }
}

unsigned int MBMBuilder::CalculateNumChemicals(const MeshState& mesh_state)
{
        // Check how many attributes represent chemicals
        // (i.e.: have the form chem_n, with n = 0 -> num_chem-1)
        vector<string> cell_attr_d = mesh_state.CellAttributeContainer().GetNames<double>();
        size_t num_chem = 0;
        vector<string> cell_chems; // Will contain the names of all attributes that are chemicals
        for_each(cell_attr_d.begin(), cell_attr_d.end(), [&](string attr_name){
                // If an attribute name starts with 'chem_' it MUST be a chemical
                if(boost::starts_with(attr_name, "chem_")) {
                        cell_chems.push_back(attr_name);
                        ++num_chem;
                }
        });

        // Check chem names collected in 'cell_chems' for consistency
        // (i.e: must be consecutive numbers from 0 to num_chem-1)
        for (size_t chem_idx = 0; chem_idx < num_chem; ++chem_idx) {
                const string chem_name = string("chem_") + boost::lexical_cast<string>(chem_idx);
                if (find(cell_chems.begin(), cell_chems.end(), chem_name) == cell_chems.end())
                        throw runtime_error("Could not find chemical: \'" + chem_name + "\'");
        }

        return num_chem;
}

void MBMBuilder::ConnectCellsToWalls(const MeshState& mesh_state)
{
        // Loop over all cells in the tissue, check which walls they contain and
        // set the appropriate wall pointers. Both cells and walls need to be constructed
        // first for this to work since access to pointers of existing walls and cells is necessary.
        auto& mesh_cells = m_mesh->GetCells();
        auto& mesh_walls = m_mesh->GetWalls();

        for (size_t cell_idx = 0; cell_idx < mesh_state.GetNumCells(); ++cell_idx) {
                const int cell_id = mesh_state.GetCellID(cell_idx);
                vector<int> cell_walls = mesh_state.GetCellWalls(cell_idx);
                // Given the wall IDs in cell_walls, find pointers to these walls and assign them to the cell
                for (const int wall_id : cell_walls) {
                        auto wall_it = find_if(mesh_walls.begin(), mesh_walls.end(), [wall_id](Wall* w) {
                                return static_cast<unsigned int>(wall_id) == w->GetIndex();
                        });
                        assert((wall_it != mesh_walls.end()) && "Cannot find the cell wall!");
                        mesh_cells[cell_id]->GetWalls().push_back(*wall_it);
                }
        }

        // Boundary polygon has walls as well
        vector<int> bp_walls = mesh_state.GetBoundaryPolygonWalls();
        // Given the wall IDs in bp_walls, find pointers to these walls and assign them to the bp
        for (const int wall_id : bp_walls) {
                auto wall_it = find_if(mesh_walls.begin(), mesh_walls.end(), [wall_id](Wall* w) {
                        return static_cast<unsigned int>(wall_id) == w->GetIndex();
                });
                assert((wall_it != mesh_walls.end()) && "Cannot find the boundary wall!");
                m_mesh->GetBoundaryPolygon()->GetWalls().push_back(*wall_it);
        }

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

}