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

#include "bio/Cell.h"
#include "bio/Mesh.h"
#include "bio/Wall.h"
#include "model/ComponentFactoryProxy.h"
#include "sim/CoreData.h"

#include <boost/property_tree/ptree.hpp>
#include <algorithm>
#include <cmath>
#include <sstream>
#include <stdexcept>

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

namespace SimPT_Sim {

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

void TransportEquations::GetDerivatives(vector<double>& derivs) const
{
        // Clear any previous entries in the vector to 0, just to be safe.
        derivs.assign(m_equation_count, 0.0);

        // Index over the concatenated vector of chem. values for cells & walls.
        size_t i = 0;

        // First part of derivs vector relates to chem. values of cells.
        for (auto c : m_mesh->GetCells()) {
                m_cell_chemistry(c, &(derivs[i]));
                i += m_effective_nchem;
        }

        // Second part of derivs vector relates to chem. values of walls.
        for (auto w : m_mesh->GetWalls()) {
                m_wall_chemistry(w, &(derivs[i]), &(derivs[i + m_effective_nchem]));

                // Transport function adds to derivatives of cell chemicals
                double* dchem_c1 = &(derivs[w->GetC1()->GetIndex() * m_effective_nchem]);
                double* dchem_c2 = &(derivs[w->GetC2()->GetIndex() * m_effective_nchem]);
                // TODO: (PK) this is not a very clean solution!
                // quick fix: dummy values to prevent end user from writing into outer space and causing a crash
                // start here if you want to implement chemical input/output into environment over boundaries
                double dummy1;
                double dummy2;
                /* Check if c1/2 is the boundary polygon */
                if (w->GetC1()->GetIndex() < 0) {
                        dchem_c1 = &dummy1;
                }
                if (w->GetC2()->GetIndex() < 0) {
                        dchem_c2 = &dummy2;
                }
                m_cell2cell_transport(w, dchem_c1, dchem_c2);
                i += 2 * m_effective_nchem;
        }
}

int TransportEquations::GetEquationCount() const
{
        return m_equation_count;
}

void TransportEquations::GetVariables(vector<double> &values) const
{
        // Index over the concatenated vector of chem. values for cells & walls.
        size_t i = 0;

        // First part of variables vector relates to chem. values of cells.
        for (auto const& cell : m_mesh->GetCells()) {
                for (size_t ch = 0; ch < static_cast<size_t>(m_effective_nchem); ch++) {
                        values[i + ch] = cell->GetChemical(ch);
                }
                i += m_effective_nchem;
        }

        // Second part of variables vector relates to chem. values of walls.
        for (auto const& wall : m_mesh->GetWalls()) {
                for (size_t ch = 0; ch < static_cast<size_t>(m_effective_nchem); ch++) {
                        values[i + ch] = wall->GetTransporters1(ch);
                }
                i += m_effective_nchem;

                for (size_t ch = 0; ch < static_cast<size_t>(m_effective_nchem); ch++) {
                        values[i + ch] = wall->GetTransporters2(ch);
                }
                i += m_effective_nchem;
        }
}

void TransportEquations::Initialize(const CoreData& cd)
{
        assert(cd.Check() && "CoreData not ok.");
        m_cd = cd;
        m_mesh = m_cd.m_mesh.get();

        const auto model_group  = m_cd.m_parameters->get<string>("model.group", "");
        const auto factory = ComponentFactoryProxy::Create(model_group);
        m_cell_chemistry = factory->CreateCellChemistry(m_cd);
        m_cell2cell_transport = factory->CreateCellToCellTransport(m_cd);
        m_wall_chemistry = factory->CreateWallChemistry(m_cd);

        const unsigned int mesh_cc = m_mesh->GetNumChemicals();
        const unsigned int model_cc = m_cd.m_parameters->get<unsigned int>("model.cell_chemical_count");
        if ((mesh_cc == 0) || (model_cc == 0) || (mesh_cc == model_cc)) {
                m_effective_nchem = min(mesh_cc, model_cc);
        } else {
                ostringstream ss;
                ss << "TransportEquations::Initialize>: mismatch in chemical counts: ";
                ss << "mesh: " << mesh_cc << " , model: " << model_cc;
                throw runtime_error(ss.str());
        }

        m_equation_count = (2 * m_mesh->GetWalls().size() + m_mesh->GetCells().size()) * m_effective_nchem;
}

void TransportEquations::SetVariables(const vector<double>& values) const
{
        // Index over the concatenated vector of chem. values for cells & walls.
        size_t i = 0;

        // First part of variables vector relates to chem. values of cells.
        for (auto cell : m_mesh->GetCells()) {
                for (size_t ch = 0; ch < static_cast<size_t>(m_effective_nchem); ch++) {
                        cell->SetChemical(ch, values[i + ch]);
                }
                i += m_effective_nchem;
        }

        // Second part of variables vector relates to chem. values of walls.
        for (const auto& wall : m_mesh->GetWalls()) {
                for (size_t ch = 0; ch < static_cast<size_t>(m_effective_nchem); ch++) {
                        wall->SetTransporters1(ch, values[i + ch]);
                }
                i += m_effective_nchem;

                for (size_t ch = 0; ch < static_cast<size_t>(m_effective_nchem); ch++) {
                        wall->SetTransporters2(ch, values[i + ch]);
                }
                i += m_effective_nchem;
        }
}

} // namespace