63 EXCEPTION(
"NagaiHondaForce is to be used with a VertexBasedCellPopulation only");
68 unsigned num_nodes = p_cell_population->
GetNumNodes();
76 bool using_target_area_modifier = p_cell_population->
Begin()->GetCellData()->HasItem(
"target area");
79 std::vector<double> element_areas(num_elements);
80 std::vector<double> element_perimeters(num_elements);
81 std::vector<double> target_areas(num_elements);
86 unsigned elem_index = elem_iter->GetIndex();
90 if (using_target_area_modifier)
92 target_areas[elem_index] = p_cell_population->
GetCellUsingLocationIndex(elem_index)->GetCellData()->GetItem(
"target area");
96 target_areas[elem_index] = mNagaiHondaTargetAreaParameter;
101 for (
unsigned node_index=0; node_index<num_nodes; node_index++)
117 c_vector<double, DIM> deformation_contribution = zero_vector<double>(DIM);
118 c_vector<double, DIM> membrane_surface_tension_contribution = zero_vector<double>(DIM);
119 c_vector<double, DIM> adhesion_contribution = zero_vector<double>(DIM);
125 for (std::set<unsigned>::iterator iter = containing_elem_indices.begin();
126 iter != containing_elem_indices.end();
131 unsigned elem_index = p_element->
GetIndex();
132 unsigned num_nodes_elem = p_element->
GetNumNodes();
139 deformation_contribution -= 2*GetNagaiHondaDeformationEnergyParameter()*(element_areas[elem_index] - target_areas[elem_index])*element_area_gradient;
142 unsigned previous_node_local_index = (num_nodes_elem+local_index-1)%num_nodes_elem;
145 unsigned next_node_local_index = (local_index+1)%num_nodes_elem;
149 double previous_edge_adhesion_parameter = GetAdhesionParameter(p_previous_node, p_this_node, *p_cell_population);
150 double next_edge_adhesion_parameter = GetAdhesionParameter(p_this_node, p_next_node, *p_cell_population);
157 adhesion_contribution -= previous_edge_adhesion_parameter*previous_edge_gradient + next_edge_adhesion_parameter*next_edge_gradient;
160 c_vector<double, DIM> element_perimeter_gradient;
161 element_perimeter_gradient = previous_edge_gradient + next_edge_gradient;
162 double cell_target_perimeter = 2*sqrt(M_PI*target_areas[elem_index]);
163 membrane_surface_tension_contribution -= 2*GetNagaiHondaMembraneSurfaceEnergyParameter()*(element_perimeters[elem_index] - cell_target_perimeter)*element_perimeter_gradient;
166 c_vector<double, DIM> force_on_node = deformation_contribution + membrane_surface_tension_contribution + adhesion_contribution;
179 std::set<unsigned> shared_elements;
180 std::set_intersection(elements_containing_nodeA.begin(),
181 elements_containing_nodeA.end(),
182 elements_containing_nodeB.begin(),
183 elements_containing_nodeB.end(),
184 std::inserter(shared_elements, shared_elements.begin()));
187 assert(!shared_elements.empty());
189 double adhesion_parameter = GetNagaiHondaCellCellAdhesionEnergyParameter();
192 if (shared_elements.size() == 1)
194 adhesion_parameter = GetNagaiHondaCellBoundaryAdhesionEnergyParameter();
197 return adhesion_parameter;