Chaste  Release::3.4
FarhadifarForce.cpp
1 /*
2 
3 Copyright (c) 2005-2016, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "FarhadifarForce.hpp"
37 
38 template<unsigned DIM>
40  : AbstractForce<DIM>(),
41  mAreaElasticityParameter(1.0), // These parameters are Case I in Farhadifar's paper
42  mPerimeterContractilityParameter(0.04),
43  mLineTensionParameter(0.12),
44  mBoundaryLineTensionParameter(0.12) // this parameter as such does not exist in Farhadifar's model.
45 {
46 }
47 
48 template<unsigned DIM>
50 {
51 }
52 
53 template<unsigned DIM>
55 {
56  // Throw an exception message if not using a VertexBasedCellPopulation
58  if (dynamic_cast<VertexBasedCellPopulation<DIM>*>(&rCellPopulation) == NULL)
59  {
60  EXCEPTION("FarhadifarForce is to be used with a VertexBasedCellPopulation only");
61  }
62 
63  // Define some helper variables
64  VertexBasedCellPopulation<DIM>* p_cell_population = static_cast<VertexBasedCellPopulation<DIM>*>(&rCellPopulation);
65  unsigned num_nodes = p_cell_population->GetNumNodes();
66  unsigned num_elements = p_cell_population->GetNumElements();
67 
68  // Begin by computing the area and perimeter of each element in the mesh, to avoid having to do this multiple times
69  std::vector<double> element_areas(num_elements);
70  std::vector<double> element_perimeters(num_elements);
71  std::vector<double> target_areas(num_elements);
72  for (typename VertexMesh<DIM,DIM>::VertexElementIterator elem_iter = p_cell_population->rGetMesh().GetElementIteratorBegin();
73  elem_iter != p_cell_population->rGetMesh().GetElementIteratorEnd();
74  ++elem_iter)
75  {
76  unsigned elem_index = elem_iter->GetIndex();
77  element_areas[elem_index] = p_cell_population->rGetMesh().GetVolumeOfElement(elem_index);
78  element_perimeters[elem_index] = p_cell_population->rGetMesh().GetSurfaceAreaOfElement(elem_index);
79  try
80  {
81  // If we haven't specified a growth modifier, there won't be any target areas in the CellData array and CellData
82  // will throw an exception that it doesn't have "target area" entries. We add this piece of code to give a more
83  // understandable message. There is a slight chance that the exception is thrown although the error is not about the
84  // target areas.
85  target_areas[elem_index] = p_cell_population->GetCellUsingLocationIndex(elem_index)->GetCellData()->GetItem("target area");
86  }
87  catch (Exception&)
88  {
89  EXCEPTION("You need to add an AbstractTargetAreaModifier to the simulation in order to use a FarhadifarForce");
90  }
91  }
92 
93  // Iterate over vertices in the cell population
94  for (unsigned node_index=0; node_index<num_nodes; node_index++)
95  {
96  Node<DIM>* p_this_node = p_cell_population->GetNode(node_index);
97 
98  /*
99  * The force on this Node is given by the gradient of the total free
100  * energy of the CellPopulation, evaluated at the position of the vertex. This
101  * free energy is the sum of the free energies of all CellPtrs in
102  * the cell population. The free energy of each CellPtr is comprised of three
103  * terms - an area deformation energy, a perimeter deformation energy
104  * and line tension energy.
105  *
106  * Note that since the movement of this Node only affects the free energy
107  * of the CellPtrs containing it, we can just consider the contributions
108  * to the free energy gradient from each of these CellPtrs.
109  */
110  c_vector<double, DIM> area_elasticity_contribution = zero_vector<double>(DIM);
111  c_vector<double, DIM> perimeter_contractility_contribution = zero_vector<double>(DIM);
112  c_vector<double, DIM> line_tension_contribution = zero_vector<double>(DIM);
113 
114  // Find the indices of the elements owned by this node
115  std::set<unsigned> containing_elem_indices = p_cell_population->GetNode(node_index)->rGetContainingElementIndices();
116 
117  // Iterate over these elements
118  for (std::set<unsigned>::iterator iter = containing_elem_indices.begin();
119  iter != containing_elem_indices.end();
120  ++iter)
121  {
122  // Get this element, its index and its number of nodes
123  VertexElement<DIM, DIM>* p_element = p_cell_population->GetElement(*iter);
124  unsigned elem_index = p_element->GetIndex();
125  unsigned num_nodes_elem = p_element->GetNumNodes();
126 
127  // Find the local index of this node in this element
128  unsigned local_index = p_element->GetNodeLocalIndex(node_index);
129 
130  // Add the force contribution from this cell's area elasticity (note the minus sign)
131  c_vector<double, DIM> element_area_gradient =
132  p_cell_population->rGetMesh().GetAreaGradientOfElementAtNode(p_element, local_index);
133  area_elasticity_contribution -= GetAreaElasticityParameter()*(element_areas[elem_index] -
134  target_areas[elem_index])*element_area_gradient;
135 
136  // Get the previous and next nodes in this element
137  unsigned previous_node_local_index = (num_nodes_elem+local_index-1)%num_nodes_elem;
138  Node<DIM>* p_previous_node = p_element->GetNode(previous_node_local_index);
139 
140  unsigned next_node_local_index = (local_index+1)%num_nodes_elem;
141  Node<DIM>* p_next_node = p_element->GetNode(next_node_local_index);
142 
143  // Compute the line tension parameter for each of these edges - be aware that this is half of the actual
144  // value for internal edges since we are looping over each of the internal edges twice
145  double previous_edge_line_tension_parameter = GetLineTensionParameter(p_previous_node, p_this_node, *p_cell_population);
146  double next_edge_line_tension_parameter = GetLineTensionParameter(p_this_node, p_next_node, *p_cell_population);
147 
148  // Compute the gradient of each these edges, computed at the present node
149  c_vector<double, DIM> previous_edge_gradient =
150  -p_cell_population->rGetMesh().GetNextEdgeGradientOfElementAtNode(p_element, previous_node_local_index);
151  c_vector<double, DIM> next_edge_gradient = p_cell_population->rGetMesh().GetNextEdgeGradientOfElementAtNode(p_element, local_index);
152 
153  // Add the force contribution from cell-cell and cell-boundary line tension (note the minus sign)
154  line_tension_contribution -= previous_edge_line_tension_parameter*previous_edge_gradient +
155  next_edge_line_tension_parameter*next_edge_gradient;
156 
157  // Add the force contribution from this cell's perimeter contractility (note the minus sign)
158  c_vector<double, DIM> element_perimeter_gradient = previous_edge_gradient + next_edge_gradient;
159  perimeter_contractility_contribution -= GetPerimeterContractilityParameter()* element_perimeters[elem_index]*
160  element_perimeter_gradient;
161  }
162 
163  c_vector<double, DIM> force_on_node = area_elasticity_contribution + perimeter_contractility_contribution + line_tension_contribution;
164  p_cell_population->GetNode(node_index)->AddAppliedForceContribution(force_on_node);
165  }
166 }
167 
168 template<unsigned DIM>
170 {
171  // Find the indices of the elements owned by each node
172  std::set<unsigned> elements_containing_nodeA = pNodeA->rGetContainingElementIndices();
173  std::set<unsigned> elements_containing_nodeB = pNodeB->rGetContainingElementIndices();
174 
175  // Find common elements
176  std::set<unsigned> shared_elements;
177  std::set_intersection(elements_containing_nodeA.begin(),
178  elements_containing_nodeA.end(),
179  elements_containing_nodeB.begin(),
180  elements_containing_nodeB.end(),
181  std::inserter(shared_elements, shared_elements.begin()));
182 
183  // Check that the nodes have a common edge
184  assert(!shared_elements.empty());
185 
186  // Since each internal edge is visited twice in the loop above, we have to use half the line tension parameter
187  // for each visit.
188  double line_tension_parameter_in_calculation = GetLineTensionParameter()/2.0;
189 
190  // If the edge corresponds to a single element, then the cell is on the boundary
191  if (shared_elements.size() == 1)
192  {
193  line_tension_parameter_in_calculation = GetBoundaryLineTensionParameter();
194  }
195 
196  return line_tension_parameter_in_calculation;
197 }
198 
199 template<unsigned DIM>
201 {
202  return mAreaElasticityParameter;
203 }
204 
205 template<unsigned DIM>
207 {
208  return mPerimeterContractilityParameter;
209 }
210 
211 template<unsigned DIM>
213 {
214  return mLineTensionParameter;
215 }
216 
217 template<unsigned DIM>
219 {
220  return mBoundaryLineTensionParameter;
221 }
222 
223 template<unsigned DIM>
224 void FarhadifarForce<DIM>::SetAreaElasticityParameter(double areaElasticityParameter)
225 {
226  mAreaElasticityParameter = areaElasticityParameter;
227 }
228 
229 template<unsigned DIM>
230 void FarhadifarForce<DIM>::SetPerimeterContractilityParameter(double perimeterContractilityParameter)
231 {
232  mPerimeterContractilityParameter = perimeterContractilityParameter;
233 }
234 
235 template<unsigned DIM>
236 void FarhadifarForce<DIM>::SetLineTensionParameter(double lineTensionParameter)
237 {
238  mLineTensionParameter = lineTensionParameter;
239 }
240 
241 template<unsigned DIM>
242 void FarhadifarForce<DIM>::SetBoundaryLineTensionParameter(double boundaryLineTensionParameter)
243 {
244  mBoundaryLineTensionParameter = boundaryLineTensionParameter;
245 }
246 
247 template<unsigned DIM>
248 void FarhadifarForce<DIM>::OutputForceParameters(out_stream& rParamsFile)
249 {
250  *rParamsFile << "\t\t\t<AreaElasticityParameter>" << mAreaElasticityParameter << "</AreaElasticityParameter>\n";
251  *rParamsFile << "\t\t\t<PerimeterContractilityParameter>" << mPerimeterContractilityParameter << "</PerimeterContractilityParameter>\n";
252  *rParamsFile << "\t\t\t<LineTensionParameter>" << mLineTensionParameter << "</LineTensionParameter>\n";
253  *rParamsFile << "\t\t\t<BoundaryLineTensionParameter>" << mBoundaryLineTensionParameter << "</BoundaryLineTensionParameter>\n";
254 
255  // Call method on direct parent class
257 }
258 
259 // Explicit instantiation
260 template class FarhadifarForce<1>;
261 template class FarhadifarForce<2>;
262 template class FarhadifarForce<3>;
263 
264 // Serialization for Boost >= 1.36
void SetAreaElasticityParameter(double areaElasticityParameter)
void OutputForceParameters(out_stream &rParamsFile)
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
Definition: Node.hpp:58
#define EXCEPTION(message)
Definition: Exception.hpp:143
double GetLineTensionParameter()
std::set< unsigned > & rGetContainingElementIndices()
Definition: Node.cpp:301
c_vector< double, SPACE_DIM > GetNextEdgeGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
MutableVertexMesh< DIM, DIM > & rGetMesh()
VertexElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
Definition: VertexMesh.hpp:668
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
Node< DIM > * GetNode(unsigned index)
unsigned GetNumNodes() const
virtual void AddForceContribution(AbstractCellPopulation< DIM > &rCellPopulation)
virtual double GetVolumeOfElement(unsigned index)
void SetPerimeterContractilityParameter(double perimeterContractilityParameter)
double GetAreaElasticityParameter()
unsigned GetIndex() const
double GetPerimeterContractilityParameter()
void SetLineTensionParameter(double lineTensionParameter)
c_vector< double, SPACE_DIM > GetAreaGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
unsigned GetNodeLocalIndex(unsigned globalIndex) const
double GetBoundaryLineTensionParameter()
VertexElementIterator GetElementIteratorEnd()
Definition: VertexMesh.hpp:675
virtual double GetSurfaceAreaOfElement(unsigned index)
virtual ~FarhadifarForce()
VertexElement< DIM, DIM > * GetElement(unsigned elementIndex)
void SetBoundaryLineTensionParameter(double boundaryLineTensionParameter)
virtual void OutputForceParameters(out_stream &rParamsFile)=0