delta_hamiltonian/PlainGC.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 "PlainGC.h"

#include "DHelper.h"
#include "bio/Cell.h"
#include "bio/GeoData.h"
#include "bio/Mesh.h"
#include "bio/NeighborNodes.h"
#include "bio/Node.h"
#include "math/Geom.h"
#include "sim/CoreData.h"
#include "util/container/circular_iterator.h"

#include <cassert>

using namespace std;
using namespace boost::property_tree;

namespace SimPT_Default {
namespace DeltaHamiltonian {

PlainGC::PlainGC(const CoreData& cd)
{
        assert( cd.Check() && "CoreData not ok in PlainGC DeltaHamiltonian");
        Initialize(cd);
}

void PlainGC::Initialize(const CoreData& cd)
{
        const auto& p = cd.m_parameters->get_child("cell_mechanics");

        m_lambda_alignment               = p.get<double>("lambda_alignment");
        m_lambda_bend                    = p.get<double>("lambda_bend");
        m_lambda_cell_length             = p.get<double>("lambda_celllength");
        m_lambda_length                  = p.get<double>("lambda_length");
        m_rp_stiffness                   = p.get<double>("relative_perimeter_stiffness");
        m_target_node_distance           = p.get<double>("target_node_distance");

}

double PlainGC::operator()(const NeighborNodes& nb, Node* n, std::array<double, 3> move)
{
        double dh = 0.0;

        const Cell* cell  = nb.GetCell();
        const auto& nb1   = nb.GetNb1();
        const auto& nb2   = nb.GetNb2();
        const array<double, 3> now_p  = *n;
        const array<double, 3> mov_p  = now_p + move;
        const array<double, 3> nb1_p  = *nb1;
        const array<double, 3> nb2_p  = *nb2;

        if (!cell->IsBoundaryPolygon()) {

                // --------------------------------------------------------------------------------------------
                // Cell area:
                // --------------------------------------------------------------------------------------------
                const auto del_area = 0.5 * ((mov_p[0] - now_p[0]) * (nb1_p[1] - nb2_p[1])
                                           + (mov_p[1] - now_p[1]) * (nb2_p[0] - nb1_p[0]));

                dh += del_area * (2.0 * (cell->GetTargetArea() - cell->GetArea()) + del_area);

                // --------------------------------------------------------------------------------------------
                // Edge length:
                // --------------------------------------------------------------------------------------------
                const bool at_b = n->IsAtBoundary();
                const auto w1   = (at_b && nb1->IsAtBoundary()) ? m_rp_stiffness : 1.0;
                const auto w2   = (at_b && nb2->IsAtBoundary()) ? m_rp_stiffness : 1.0;
                const auto old1 = Norm(now_p - nb1_p) - m_target_node_distance;
                const auto old2 = Norm(now_p - nb2_p) - m_target_node_distance;
                const auto new1 = Norm(mov_p - nb1_p) - m_target_node_distance;
                const auto new2 = Norm(mov_p - nb2_p) - m_target_node_distance;

                dh += m_lambda_length * (w1 * (new1 * new1 - old1 * old1) + w2 * (new2 * new2 - old2 * old2));

                // --------------------------------------------------------------------------------------------
                // Cell length and alignment:
                // --------------------------------------------------------------------------------------------
                if ((abs(m_lambda_cell_length) > 1.0e-7) || (abs(m_lambda_alignment) > 1.0e-7)) {
                        const auto old_geo = cell->GetGeoData();
                        const auto new_geo = ModifyForMove(old_geo, nb1_p, nb2_p, now_p, mov_p);
                        const auto old_tup = old_geo.GetEllipseAxes();
                        const auto new_tup = new_geo.GetEllipseAxes();

                        // ------------------------------------------------------------------------------------
                        // CellLength:
                        // ------------------------------------------------------------------------------------
                        if (abs(m_lambda_cell_length) > 1.0e-7) {
                                const double t_old = cell->GetTargetLength() - get<0>(old_tup);
                                const double t_new = cell->GetTargetLength() - get<0>(new_tup);

                                dh += m_lambda_cell_length * (t_new * t_new - t_old * t_old);
                        }

                        // ------------------------------------------------------------------------------------
                        // Alignment:
                        // ------------------------------------------------------------------------------------
                        if (abs(m_lambda_alignment) > 1.0e-7) {
                                dh += m_lambda_alignment * DHelper::AlignmentTerm(get<1>(old_tup), get<1>(new_tup));
                        }
                }
        }

        // --------------------------------------------------------------------------------------------
        // Vertex bending:
        // --------------------------------------------------------------------------------------------
        if (abs(m_lambda_bend) > 1.0e-7) {
                dh += m_lambda_bend * DHelper::BendingTerm(nb1_p, nb2_p, now_p, mov_p);
        }

        return dh;
}

} // namespace
} // namespace