DHelper.h

Go to the documentation of this file.

#ifndef SIMPT_DEFAULT_DELTA_HAMILTONIAN_DHELPER_H_INCLUDED
#define SIMPT_DEFAULT_DELTA_HAMILTONIAN_DHELPER_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 "bio/Wall.h"
#include "math/array3.h"
#include "math/Geom.h"
#include "util/container/circular_iterator.h"

#include <array>

using namespace std;
using namespace boost::property_tree;

namespace SimPT_Default {
namespace DeltaHamiltonian {

using namespace SimPT_Sim;

class DHelper
{
public:
        static double AlignmentTerm(
                const std::array<double, 3>& old_axis, const std::array<double, 3>& new_axis)
        {
                // orth_alignment orthogonal to intended alignment.
                const auto old_a = Normalize(old_axis);
                const auto new_a = Normalize(new_axis);
                const array<double, 3> orth_alignment {{-1.0, 0.0, 0.0}};
                const double alignment_old  = InnerProduct(old_a, orth_alignment);
                const double alignment_new  = InnerProduct(new_a, orth_alignment);

                return (alignment_old - alignment_new);
        }

        static double BendingTerm(
                const array<double, 3>& nb1_p, const array<double, 3>& nb2_p,
                const array<double, 3>& now_p, const array<double, 3>& mov_p)
        {
                auto          t1  = Geom::CircumCircle(nb1_p[0], nb1_p[1], now_p[0], now_p[1], nb2_p[0], nb2_p[1]);
                double const  r1  = get<2>(t1);
                auto          t2  = Geom::CircumCircle(nb1_p[0], nb1_p[1], mov_p[0], mov_p[1], nb2_p[0], nb2_p[1]);
                double const  r2  = get<2>(t2);
                assert(r1 > 0.0 && r2 > 0.0 && "Circumcenter radius singular!");

                return pow((1.0 / r2 - 1.0 / r1), 2);
        }

        static double ElasticWallTerm(
                const Cell* cell, Node* n, const std::list<Wall*>& node_owning_walls, double stiffness,
                const std::array<double, 3>& nb1_p, const std::array<double, 3>& nb2_p,
                const std::array<double, 3>& now_p, const std::array<double, 3>& mov_p)
        {
                double contrib = 0.0;
                for (auto& wall : node_owning_walls) {

                        const auto& c1 = wall->GetC1();
                        const auto& c2 = wall->GetC2();
                        if (c1 == cell || c2 == cell) {

                                const bool check    = (cell == c1 || cell->IsBoundaryPolygon());
                                const auto wall_beg = check ? wall->GetN1() : wall->GetN2();
                                const auto wall_end = check ? wall->GetN2() : wall->GetN1();
                                const auto stiff    = (c1->IsBoundaryPolygon() || c2->IsBoundaryPolygon()) ? stiffness * wall->GetStrength() : 1.0 * wall->GetStrength();
                                const auto rest_l   = wall->GetRestLength();
                                const auto old_l    = wall->GetLength();

                                double new_l = old_l;
                                new_l += (n != wall_beg) ? Norm(mov_p - nb1_p) - Norm(now_p - nb1_p) : 0.0;
                                new_l += (n != wall_end) ? Norm(nb2_p - mov_p) - Norm(nb2_p - now_p) : 0.0;
                                contrib += stiff * (pow(new_l - rest_l, 2) - pow(old_l - rest_l, 2)) / (2.0 * rest_l);
                        }
                }
                return contrib;
        }
};

} // namespace
} // namespace

#endif // end-of-include-guard