delta_hamiltonian/ElasticWall.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 "DHelper.h"
#include "ElasticWall.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 {

using namespace SimPT_Sim;

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

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

        m_elastic_modulus                = p.get<double>("elastic_modulus");
        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", 0.0);
        m_rp_stiffness                   = p.get<double>("relative_perimeter_stiffness");
        m_target_node_distance           = p.get<double>("target_node_distance");
}

double ElasticWall::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);

                // --------------------------------------------------------------------------------------------
                // Elastic wall:
                // --------------------------------------------------------------------------------------------
                dh += m_elastic_modulus * DHelper::ElasticWallTerm(
                        cell, n, m_mesh->GetNodeOwningWalls(n), m_rp_stiffness, nb1_p, nb2_p, now_p, mov_p);

                // --------------------------------------------------------------------------------------------
                // Edge length:
                // --------------------------------------------------------------------------------------------
                if (abs(m_lambda_length) > 1.0e-7) {
                        const bool at_b   = n->IsAtBoundary();
                        const double w1   = (at_b && nb1->IsAtBoundary()) ? m_rp_stiffness : 1.0;
                        const double w2   = (at_b && nb2->IsAtBoundary()) ? m_rp_stiffness : 1.0;
                        const double old1 = Norm(now_p - nb1_p) - m_target_node_distance;
                        const double old2 = Norm(now_p - nb2_p) - m_target_node_distance;
                        const double new1 = Norm(mov_p - nb1_p) - m_target_node_distance;
                        const double 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 (note: new_orth_axis orthogonalized through definition of its components)
                        // ------------------------------------------------------------------------------------
                        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