50 EXCEPTION(
"RepulsionForce is to be used with a NodeBasedCellPopulation only");
55 for (
typename std::vector< std::pair<
Node<DIM>*,
Node<DIM>* > >::iterator iter = r_node_pairs.begin();
56 iter != r_node_pairs.end();
59 std::pair<Node<DIM>*,
Node<DIM>* > pair = *iter;
65 const c_vector<double, DIM>& r_node_a_location = p_node_a->
rGetLocation();
66 const c_vector<double, DIM>& r_node_b_location = p_node_b->
rGetLocation();
69 double node_a_radius = p_node_a->
GetRadius();
70 double node_b_radius = p_node_b->
GetRadius();
73 c_vector<double, DIM> unit_difference;
75 unit_difference = (
static_cast<NodeBasedCellPopulation<DIM>*
>(&rCellPopulation))->rGetMesh().GetVectorFromAtoB(r_node_a_location, r_node_b_location);
78 double rest_length = node_a_radius+node_b_radius;
80 if (norm_2(unit_difference) < rest_length)
83 c_vector<double, DIM> force = this->CalculateForceBetweenNodes(p_node_a->
GetIndex(), p_node_b->
GetIndex(), rCellPopulation);
84 c_vector<double, DIM> negative_force = -1.0 * force;
85 for (
unsigned j=0; j<DIM; j++)
87 assert(!std::isnan(force[j]));