66 c_vector<double, DIM> unit_vector;
70 cell_iter != rCellPopulation.
End();
79 const c_vector<double, DIM>& r_node_i_location = p_node_i->
rGetLocation();
82 double radius_of_cell_i = p_node_i->
GetRadius();
84 double delta_V_c = 0.0;
85 c_vector<double, DIM> dVAdd_vector = zero_vector<double>(DIM);
91 for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
92 iter != neighbouring_node_indices.end();
98 const c_vector<double, DIM>& r_node_j_location = p_node_j->
rGetLocation();
101 unit_vector = r_node_j_location - r_node_i_location;
104 double dij = norm_2(unit_vector);
109 double radius_of_cell_j = p_node_j->
GetRadius();
112 if (dij < radius_of_cell_i + radius_of_cell_j)
115 double xij = 0.5*(radius_of_cell_i*radius_of_cell_i - radius_of_cell_j*radius_of_cell_j + dij*dij)/dij;
116 double dxijdd = 1.0 - xij/dij;
117 double dVAdd = M_PI*dxijdd*(5.0*pow(radius_of_cell_i,2.0) + 3.0*pow(xij,2.0) - 8.0*radius_of_cell_i*xij)/3.0;
119 dVAdd_vector += dVAdd*unit_vector;
122 delta_V_c += M_PI*pow(radius_of_cell_i - xij,2.0)*(2*radius_of_cell_i - xij)/3.0;
126 double V_A = 4.0/3.0*M_PI*pow(radius_of_cell_i,3.0) - delta_V_c;
136 c_vector<double, DIM> applied_force = -mCompressionEnergyParameter/V_T*(V_T - V_A)*dVAdd_vector;