50 std::vector< std::vector<double> > extremal_normal_forces;
51 std::vector<double> minimum_normal_forces(num_nodes);
52 std::vector<double> maximum_normal_forces(num_nodes);
54 for (
unsigned i=0; i<num_nodes; i++)
59 double minimum_normal_force_for_node_i = DBL_MAX;
60 double maximum_normal_force_for_node_i = -DBL_MAX;
62 for (
unsigned j=0; j<extremal_angles.size(); j++)
66 if (current_normal_force > maximum_normal_force_for_node_i)
68 maximum_normal_force_for_node_i = current_normal_force;
70 if (current_normal_force < minimum_normal_force_for_node_i)
72 minimum_normal_force_for_node_i = current_normal_force;
76 assert( minimum_normal_force_for_node_i <= maximum_normal_force_for_node_i);
78 minimum_normal_forces[i] = minimum_normal_force_for_node_i;
79 maximum_normal_forces[i] = maximum_normal_force_for_node_i;
82 extremal_normal_forces.push_back(minimum_normal_forces);
83 extremal_normal_forces.push_back(maximum_normal_forces);
85 return extremal_normal_forces;
91 std::ostringstream time_string;
93 std::string results_directory = simulationOutputDirectory +
"/results_from_time_" + time_string.str();
98 (*mpVizStressResultsFile) << time <<
"\t";
112 global_index = (
double) i;
114 x = r_mesh.
GetNode(i)->rGetLocation()[0];
115 y = r_mesh.
GetNode(i)->rGetLocation()[1];
117 minimum = extremal_normal_forces[0][i];
118 maximum = extremal_normal_forces[1][i];
120 (*mpVizStressResultsFile) << global_index <<
" " << x <<
" " << y <<
" " << minimum <<
" " << maximum <<
" ";
123 (*mpVizStressResultsFile) <<
"\n";
133 double tangential_force = 0.0;
134 double normal_force = 0.0;
137 c_vector<double,2> unit_vec_between_nodes(2);
139 for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
140 iter != neighbouring_node_indices.end();
146 assert(alpha <= M_PI);
147 assert(alpha > -M_PI);
149 if (sin(alpha-theta) > DBL_EPSILON)
152 c_vector<double,2> force_between_nodes = zero_vector<double>(2);
159 force_between_nodes += (*force_iter)->CalculateForceBetweenNodes(index, *iter,
mrCellPopulation);
162 unit_vec_between_nodes[0] = cos(alpha);
163 unit_vec_between_nodes[1] = sin(alpha);
165 double plusminus_norm_force = inner_prod(force_between_nodes,unit_vec_between_nodes);
166 tangential_force += plusminus_norm_force * cos(alpha-theta);
167 normal_force += plusminus_norm_force * sin(alpha-theta);
171 std::vector<double> ret(2);
172 ret[0] = tangential_force;
173 ret[1] = normal_force;
184 std::vector<double> sampling_angles(4*neighbouring_node_indices.size());
188 for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
189 iter != neighbouring_node_indices.end();
195 double alpha_minus_epsilon = alpha -
mEpsilon;
196 double alpha_plus_epsilon = alpha +
mEpsilon;
197 double alpha_plus_pi_minus_epsilon = alpha + M_PI -
mEpsilon;
198 double alpha_plus_pi_plus_epsilon = alpha + M_PI +
mEpsilon;
201 if (alpha_minus_epsilon <= -M_PI)
203 alpha_minus_epsilon += 2*M_PI;
205 sampling_angles[i] = alpha_minus_epsilon;
207 assert(sampling_angles[i] <= M_PI);
208 assert(sampling_angles[i] > -M_PI);
211 if (alpha_plus_epsilon > M_PI)
213 alpha_plus_epsilon -= 2*M_PI;
215 sampling_angles[i] = alpha_plus_epsilon;
217 assert(sampling_angles[i] <= M_PI);
218 assert(sampling_angles[i] > -M_PI);
221 if (alpha_plus_pi_minus_epsilon > M_PI)
223 alpha_plus_pi_minus_epsilon -= 2*M_PI;
225 sampling_angles[i] = alpha_plus_pi_minus_epsilon;
227 assert(sampling_angles[i] <= M_PI);
228 assert(sampling_angles[i] > -M_PI);
231 if (alpha_plus_pi_plus_epsilon > M_PI)
233 alpha_plus_pi_plus_epsilon -= 2*M_PI;
235 sampling_angles[i] = alpha_plus_pi_plus_epsilon;
237 assert(sampling_angles[i] <= M_PI);
238 assert(sampling_angles[i] > -M_PI);
242 sort(sampling_angles.begin(), sampling_angles.end());
243 return sampling_angles;
250 assert(angle1 < angle2);
252 double tolerance = 1e-5;
253 double previous_angle;
254 double current_error;
255 double current_angle = angle1;
257 current_error = angle2 - angle1;
258 std::vector<double> current_ft_and_fn(2);
260 while (current_error > tolerance)
262 previous_angle = current_angle;
264 current_angle -= current_ft_and_fn[0]/current_ft_and_fn[1];
265 current_error = fabs(current_angle - previous_angle);
268 assert(current_angle>angle1 && current_angle<angle2);
269 assert(current_error < tolerance);
271 return current_angle;
276 std::vector<double> extremal_angles;
277 std::vector<double> ft_and_fn(2);
278 std::vector<double> tangential_force(samplingAngles.size());
280 for (
unsigned i=0; i<samplingAngles.size(); i++)
283 tangential_force[i] = ft_and_fn[0];
286 unsigned n = samplingAngles.size()-1;
288 for (
unsigned i=0; i<n; i++)
290 if ((tangential_force[i%n]>0 && tangential_force[(i+1)%n]<0)
291 || (tangential_force[i%n]<0 && tangential_force[(i+1)%n]>0))
293 double next_extremal_angle;
299 samplingAngles[i%n] -= 2*M_PI;
302 if (samplingAngles[(i+1)%n] - samplingAngles[i%n] < 2*
mEpsilon + 1e-6 )
306 next_extremal_angle = 0.5*(samplingAngles[(i+1)%n] + samplingAngles[i%n]);
311 next_extremal_angle =
GetLocalExtremum(index, samplingAngles[i%n], samplingAngles[(i+1)%n]);
314 if (next_extremal_angle <= -M_PI)
316 next_extremal_angle += 2*M_PI;
318 assert(next_extremal_angle>-M_PI && next_extremal_angle<=M_PI);
319 extremal_angles.push_back(next_extremal_angle);
323 return extremal_angles;