delta_hamiltonian/ModifiedGC.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 "ModifiedGC.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 {

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

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

        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 ModifiedGC::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;
        const auto& owningWalls = m_mesh->GetNodeOwningWalls(n);

        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]));
                const auto old_area = cell->GetArea();
                const auto tar_area = cell->GetTargetArea();
                const auto old_fact = (old_area - tar_area) / old_area;
                const auto new_fact = (old_area - del_area - tar_area) / (old_area - del_area);

                dh += 1000 * (new_fact * new_fact - old_fact * old_fact);

                // --------------------------------------------------------------------------------------------
                // Wall length:
                // --------------------------------------------------------------------------------------------
                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;

                const bool at_b = n->IsAtBoundary();
                double w1 = (at_b && nb1->IsAtBoundary()) ? m_rp_stiffness : 1.0;
                double w2 = (at_b && nb2->IsAtBoundary()) ? m_rp_stiffness : 1.0;

                for (auto& wall : owningWalls) {
                        if (wall->GetC1() == cell || wall->GetC2() == cell) {
                                const auto w  = wall->GetStrength();
                                const auto n1 = wall->GetN1();
                                const auto n2 = wall->GetN2();

                                if ( ((n1 == n) && (n2 == nb1)) || ((n2 == n) && (n1 == nb1)) ) {
                                        w1 = (at_b && nb1->IsAtBoundary()) ? m_rp_stiffness * w : w;
                                }
                                if ( ((n1 == n) && (n2 == nb2)) || ((n2 == n) && (n1 == nb2)) ) {
                                        w2 = (at_b && nb2->IsAtBoundary()) ? m_rp_stiffness * w : w;
                                }
                        }
                }

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

        return dh;
}

} // namespace
} // namespace