55 if constexpr (DIM == 2)
60 EXCEPTION(
"WelikyOsterForce is to be used with a VertexBasedCellPopulation only");
76 unsigned element_index = element_iter->GetIndex();
83 double deformation_coefficient = GetWelikyOsterAreaParameter()/element_area;
92 double membrane_surface_tension_coefficient = GetWelikyOsterPerimeterParameter()*element_perimeter;
96 unsigned num_nodes = element_iter->GetNumNodes();
97 for (
unsigned node_local_index = 0; node_local_index < num_nodes; node_local_index++)
99 unsigned node_global_index = element_iter->GetNodeGlobalIndex(node_local_index);
101 c_vector<double, DIM> current_node = element_iter->GetNodeLocation(node_local_index);
102 c_vector<double, DIM> next_node = element_iter->GetNodeLocation((node_local_index + 1)%(element_iter->GetNumNodes()));
103 c_vector<double, DIM> previous_node = element_iter->GetNodeLocation((node_local_index + element_iter->GetNumNodes() - 1)%(element_iter->GetNumNodes()));
105 c_vector<double, DIM> clockwise_unit_vector = p_cell_population->
rGetMesh().
GetVectorFromAtoB(current_node, previous_node);
106 clockwise_unit_vector /= norm_2(clockwise_unit_vector);
108 c_vector<double, DIM> anti_clockwise_unit_vector = p_cell_population->
rGetMesh().
GetVectorFromAtoB(current_node, next_node);
109 anti_clockwise_unit_vector /= norm_2(anti_clockwise_unit_vector);
112 c_vector<double, DIM> outward_normal = -0.5*clockwise_unit_vector - 0.5*anti_clockwise_unit_vector;
113 outward_normal /= norm_2(outward_normal);
115 c_vector<double, DIM> deformation_contribution = deformation_coefficient * outward_normal;
117 c_vector<double, DIM> membrane_surface_tension_contribution = membrane_surface_tension_coefficient * (clockwise_unit_vector + anti_clockwise_unit_vector);
119 c_vector<double, DIM> force_on_node = deformation_contribution + membrane_surface_tension_contribution;