VPTissue Reference Manual
delta_hamiltonian/ElasticWall.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2011-2016 Universiteit Antwerpen
3  *
4  * Licensed under the EUPL, Version 1.1 or as soon they will be approved by
5  * the European Commission - subsequent versions of the EUPL (the "Licence");
6  * You may not use this work except in compliance with the Licence.
7  * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
8  *
9  * Unless required by applicable law or agreed to in writing, software
10  * distributed under the Licence is distributed on an "AS IS" basis,
11  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12  * See the Licence for the specific language governing
13  * permissions and limitations under the Licence.
14  */
20 #include "DHelper.h"
21 #include "ElasticWall.h"
22 
23 #include "bio/Cell.h"
24 #include "bio/GeoData.h"
25 #include "bio/Mesh.h"
26 #include "bio/NeighborNodes.h"
27 #include "bio/Node.h"
28 #include "math/Geom.h"
29 #include "sim/CoreData.h"
31 
32 #include <cassert>
33 
34 using namespace std;
35 using namespace boost::property_tree;
36 
37 namespace SimPT_Default {
38 namespace DeltaHamiltonian {
39 
40 using namespace SimPT_Sim;
41 
42 ElasticWall::ElasticWall(const CoreData& cd)
43 {
44  assert( cd.Check() && "CoreData not ok in ElasticWall DeltaHamiltonian");
45  Initialize(cd);
46 }
47 
48 void ElasticWall::Initialize(const CoreData& cd)
49 {
50  m_mesh = cd.m_mesh.get();
51  const auto& p = cd.m_parameters->get_child("cell_mechanics");
52 
53  m_elastic_modulus = p.get<double>("elastic_modulus");
54  m_lambda_alignment = p.get<double>("lambda_alignment");
55  m_lambda_bend = p.get<double>("lambda_bend");
56  m_lambda_cell_length = p.get<double>("lambda_celllength");
57  m_lambda_length = p.get<double>("lambda_length", 0.0);
58  m_rp_stiffness = p.get<double>("relative_perimeter_stiffness");
59  m_target_node_distance = p.get<double>("target_node_distance");
60 }
61 
62 double ElasticWall::operator()(const NeighborNodes& nb, Node* n, std::array<double, 3> move)
63 {
64  double dh = 0.0;
65 
66  const Cell* cell = nb.GetCell();
67  const auto& nb1 = nb.GetNb1();
68  const auto& nb2 = nb.GetNb2();
69  const array<double, 3> now_p = *n;
70  const array<double, 3> mov_p = now_p + move;
71  const array<double, 3> nb1_p = *nb1;
72  const array<double, 3> nb2_p = *nb2;
73 
74  if (!cell->IsBoundaryPolygon()) {
75 
76  // --------------------------------------------------------------------------------------------
77  // Cell area:
78  // --------------------------------------------------------------------------------------------
79  const auto del_area = 0.5 * ((mov_p[0] - now_p[0]) * (nb1_p[1] - nb2_p[1])
80  + (mov_p[1] - now_p[1]) * (nb2_p[0] - nb1_p[0]));
81 
82  dh += del_area * (2.0 * (cell->GetTargetArea() - cell->GetArea()) + del_area);
83 
84  // --------------------------------------------------------------------------------------------
85  // Elastic wall:
86  // --------------------------------------------------------------------------------------------
87  dh += m_elastic_modulus * DHelper::ElasticWallTerm(
88  cell, n, m_mesh->GetNodeOwningWalls(n), m_rp_stiffness, nb1_p, nb2_p, now_p, mov_p);
89 
90  // --------------------------------------------------------------------------------------------
91  // Edge length:
92  // --------------------------------------------------------------------------------------------
93  if (abs(m_lambda_length) > 1.0e-7) {
94  const bool at_b = n->IsAtBoundary();
95  const double w1 = (at_b && nb1->IsAtBoundary()) ? m_rp_stiffness : 1.0;
96  const double w2 = (at_b && nb2->IsAtBoundary()) ? m_rp_stiffness : 1.0;
97  const double old1 = Norm(now_p - nb1_p) - m_target_node_distance;
98  const double old2 = Norm(now_p - nb2_p) - m_target_node_distance;
99  const double new1 = Norm(mov_p - nb1_p) - m_target_node_distance;
100  const double new2 = Norm(mov_p - nb2_p) - m_target_node_distance;
101 
102  dh += m_lambda_length * (w1 * (new1 * new1 - old1 * old1) + w2 * (new2 * new2 - old2 * old2));
103  }
104 
105  // --------------------------------------------------------------------------------------------
106  // Cell length and alignment:
107  // --------------------------------------------------------------------------------------------
108  if ((abs(m_lambda_cell_length) > 1.0e-7) || (abs(m_lambda_alignment) > 1.0e-7)) {
109  const auto old_geo = cell->GetGeoData();
110  const auto new_geo = ModifyForMove(old_geo, nb1_p, nb2_p, now_p, mov_p);
111  const auto old_tup = old_geo.GetEllipseAxes();
112  const auto new_tup = new_geo.GetEllipseAxes();
113 
114  // ------------------------------------------------------------------------------------
115  // CellLength:
116  // ------------------------------------------------------------------------------------
117  if (abs(m_lambda_cell_length) > 1.0e-7) {
118  const double t_old = cell->GetTargetLength() - get<0>(old_tup);
119  const double t_new = cell->GetTargetLength() - get<0>(new_tup);
120 
121  dh += m_lambda_cell_length * (t_new * t_new - t_old * t_old);
122  }
123 
124  // ------------------------------------------------------------------------------------
125  // Alignment (note: new_orth_axis orthogonalized through definition of its components)
126  // ------------------------------------------------------------------------------------
127  if (abs(m_lambda_alignment) > 1.0e-7) {
128  dh += m_lambda_alignment * DHelper::AlignmentTerm(get<1>(old_tup), get<1>(new_tup));
129  }
130  }
131  }
132 
133  // --------------------------------------------------------------------------------------------
134  // Vertex bending:
135  // --------------------------------------------------------------------------------------------
136  if (abs(m_lambda_bend) > 1.0e-7) {
137  dh += m_lambda_bend * DHelper::BendingTerm(nb1_p, nb2_p, now_p, mov_p);
138  }
139 
140  return dh;
141 }
142 
143 } // namespace
144 } // namespace
Core data with mesh, parameters, random engine and time data.
Definition: CoreData.h:38
STL namespace.
A cell contains walls and nodes.
Definition: Cell.h:48
Node in cell wall.
Definition: Node.h:39
Cell * GetCell() const
Return the cell of this Neighbor pair.
Core data used during model execution.
Combo header for circular iterator.
Interface for NeighborNodes.
Interface/Implementation for GeoData.
Namespace for components of the Default model group.
Namespace for the core simulator.
Interface for Cell.
Node * GetNb1() const
Return first node of this Neighbor pair.
Structure of neighboring nodes: two neighboring nodes from standpoint of a given cell with an orderin...
Definition: NeighborNodes.h:34
Interface for Node.
Node * GetNb2() const
Return second node of this Neighbor pair.
Interface for Geom.
Interface/Implementation for DeltaHamiltonian::DHelper.
GeoData GetGeoData() const
Return GeData (area, centroid, area moment of inertia).
Definition: Cell.cpp:225
Interface for DeltaHamiltonian::ElasticWall.
double GetArea() const
Return the area of the cell.
Definition: Cell.cpp:178
bool Check() const
Verify all pointers non-null.
Definition: CoreData.h:53
Interface for Mesh.