60 unsigned system_size = this->GetNumNodes()*DIM;
75 cell_iter != this->End();
79 unsigned global_node_index = this->GetLocationIndexUsingCell((*cell_iter));
82 unsigned node_index = this->rGetMesh().SolveNodeMapping(global_node_index);
84 Node<DIM>* p_node_i = this->GetNode(global_node_index);
87 const c_vector<double, DIM>& r_node_i_location = p_node_i->
rGetLocation();
90 double radius_of_cell_i = p_node_i->
GetRadius();
93 double damping_const = this->GetDampingConstant(global_node_index);
98 std::set<unsigned> neighbouring_node_indices = this->GetNeighbouringNodeIndices(global_node_index);
100 for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
101 iter != neighbouring_node_indices.end();
104 unsigned neighbour_node_global_index = *iter;
106 unsigned neighbour_node_index = this->rGetMesh().SolveNodeMapping(neighbour_node_global_index);
111 Node<DIM>* p_node_j = this->GetNode(neighbour_node_global_index);
114 const c_vector<double, DIM>& r_node_j_location = p_node_j->
rGetLocation();
117 c_vector<double, DIM> unit_vector = r_node_j_location - r_node_i_location;
120 double dij = norm_2(unit_vector);
125 double radius_of_cell_j = p_node_j->
GetRadius();
127 if (dij < radius_of_cell_i + radius_of_cell_j)
130 double xij = 0.5*(radius_of_cell_i*radius_of_cell_i - radius_of_cell_j*radius_of_cell_j + dij*dij)/dij;
132 Aij = M_PI*(radius_of_cell_i*radius_of_cell_i - xij*xij);
135 for (
unsigned i=0; i<DIM; i++)
144 for (
unsigned i=0; i<DIM; i++)
152 c_vector<double, DIM> current_location;
153 c_vector<double, DIM> forces;
154 current_location = this->GetNode(global_node_index)->rGetLocation();
155 forces = this->GetNode(global_node_index)->rGetAppliedForce();
157 for (
unsigned i=0; i<DIM; i++)
174 cell_iter != this->End();
178 unsigned global_node_index = this->GetLocationIndexUsingCell((*cell_iter));
180 unsigned node_index = this->rGetMesh().SolveNodeMapping(global_node_index);
182 c_vector<double, DIM> new_node_location;
185 for (
unsigned i=0; i<DIM; i++)
187 new_node_location(i) = soln_next_timestep_repl[DIM*node_index+i];
194 this->SetNode(global_node_index, new_point);