HHelper.h

Go to the documentation of this file.

#ifndef HAMILTONIAN_HHELPER_H_INCLUDED
#define HAMILTONIAN_HHELPER_H_INCLUDED
/*
 * 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 "bio/Cell.h"
#include "bio/Node.h"
#include "math/Geom.h"
#include "util/container/circular_iterator.h"

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

namespace SimPT_Default {
namespace Hamiltonian {

class HHelper
{
public:
        static double BendingTerm(Cell* cell)
        {
                // TODO: 1.0 / 10.0 hidden parameter: curvature of decagon with edge length of 6.18034.
                double contrib = 0.0;

                const auto cit_start = cell->GetNodes().begin();
                auto cit = make_const_circular(cell->GetNodes());
                do {
                        array<double, 3>  ncurr = *(*cit);
                        array<double, 3>  nprev = *(*prev(cit));
                        array<double, 3>  nnext = *(*next(cit));
                        const auto t = Geom::CircumCircle(nprev[0], nprev[1], ncurr[0], ncurr[1], nnext[0], nnext[1]);
                        const auto r = get<2>(t);
                        contrib += pow((1.0 / r - 1.0 / 10.0), 2);
                } while (++cit != cit_start);

                return contrib;
        }

        static double CellLengthTerm(Cell* cell)
        {
                return pow(get<0>(cell->GetGeoData().GetEllipseAxes()) - cell->GetTargetLength(), 2);
        }

        static double EdgeTerm(Cell* cell, double stiffness, double distance)
        {
                double  contrib      = 0.0;
                const auto cit_start = cell->GetNodes().begin();

                auto cit = make_const_circular(cell->GetNodes());
                do {
                        const auto ncurr  = *cit;
                        const auto nnext  = *next(cit);
                        const auto weight = ncurr->IsAtBoundary() && nnext->IsAtBoundary() ? stiffness : 1.0;
                        contrib += weight * pow(distance - Norm(*ncurr - *nnext), 2);
                } while (++cit != cit_start);

                return contrib;
        }

        static double ElasticWallTerm(Cell* cell, double stiffness)
        {
                double  contrib = 0.0;
                const auto wit_start = cell->GetWalls().begin();

                auto cit = make_const_circular(cell->GetWalls());
                do {
                        const auto& wall = *(*cit);
                        const double w_stiff
                                = (wall.GetC1()->IsBoundaryPolygon() || wall.GetC2()->IsBoundaryPolygon()) ? stiffness : 1.0;
                        contrib += w_stiff
                                * pow(wall.GetLength() - wall.GetRestLength(), 2) / (2. * wall.GetRestLength());
                } while (++cit != wit_start);

                return contrib;
        }

        static double WallLengthTerm(Cell* cell, double stiffness, double distance)
        {
                double contrib = 0.0;
                const auto wit_start = cell->GetWalls().begin();

                auto wit = make_const_circular(cell->GetWalls());
                do {
                        const auto& wall = *(*wit);
                        const auto s = wall.IsAtBoundary() ? stiffness : 1.0;
                        contrib += s * wall.GetStrength() * pow(distance - wall.GetLength(), 2);
                } while (++wit != wit_start);

                return contrib;
        }
};

} // namespace
} // namespace

#endif // end-of-include-guard