Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
DiscreteSystemForceCalculator.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "DiscreteSystemForceCalculator.hpp"
37
39 std::vector<boost::shared_ptr<AbstractTwoBodyInteractionForce<2> > > forceCollection)
40 : mrCellPopulation(rCellPopulation),
41 mForceCollection(forceCollection),
42 mEpsilon(0.01)
43{
44}
45
47{
48 unsigned num_nodes = mrCellPopulation.GetNumNodes();
49
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);
53
54 for (unsigned i=0; i<num_nodes; i++)
55 {
56 std::vector<double> sampling_angles = GetSamplingAngles(i);
57 std::vector<double> extremal_angles = GetExtremalAngles(i, sampling_angles);
58
59 double minimum_normal_force_for_node_i = DBL_MAX;
60 double maximum_normal_force_for_node_i = -DBL_MAX;
61
62 for (unsigned j=0; j<extremal_angles.size(); j++)
63 {
64 double current_normal_force = CalculateFtAndFn(i, extremal_angles[j])[1];
65
66 if (current_normal_force > maximum_normal_force_for_node_i)
67 {
68 maximum_normal_force_for_node_i = current_normal_force;
69 }
70 if (current_normal_force < minimum_normal_force_for_node_i)
71 {
72 minimum_normal_force_for_node_i = current_normal_force;
73 }
74 }
75
76 assert( minimum_normal_force_for_node_i <= maximum_normal_force_for_node_i);
77
78 minimum_normal_forces[i] = minimum_normal_force_for_node_i;
79 maximum_normal_forces[i] = maximum_normal_force_for_node_i;
80 }
81
82 extremal_normal_forces.push_back(minimum_normal_forces);
83 extremal_normal_forces.push_back(maximum_normal_forces);
84
85 return extremal_normal_forces;
86}
87
88void DiscreteSystemForceCalculator::WriteResultsToFile(std::string simulationOutputDirectory)
89{
90 double time = SimulationTime::Instance()->GetTime();
91 std::ostringstream time_string;
92 time_string << time;
93 std::string results_directory = simulationOutputDirectory + "/results_from_time_" + time_string.str();
94
95 OutputFileHandler output_file_handler2(results_directory+"/", false);
96 mpVizStressResultsFile = output_file_handler2.OpenOutputFile("results.vizstress");
97
98 (*mpVizStressResultsFile) << time << "\t";
99
100 double global_index;
101 double x;
102 double y;
103 double minimum;
104 double maximum;
105
107
108 std::vector< std::vector<double> > extremal_normal_forces = CalculateExtremalNormalForces();
109
110 for (unsigned i=0; i<r_mesh.GetNumNodes(); i++)
111 {
112 global_index = (double) i;
113
114 x = r_mesh.GetNode(i)->rGetLocation()[0];
115 y = r_mesh.GetNode(i)->rGetLocation()[1];
116
117 minimum = extremal_normal_forces[0][i];
118 maximum = extremal_normal_forces[1][i];
119
120 (*mpVizStressResultsFile) << global_index << " " << x << " " << y << " " << minimum << " " << maximum << " ";
121 }
122
123 (*mpVizStressResultsFile) << "\n";
124 mpVizStressResultsFile->close();
125}
126
127std::vector<double> DiscreteSystemForceCalculator::CalculateFtAndFn(unsigned index, double theta)
128{
130
131 std::set<unsigned> neighbouring_node_indices = mrCellPopulation.GetNeighbouringNodeIndices(index);
132
133 double tangential_force = 0.0;
134 double normal_force = 0.0;
135 double alpha;
136
137 c_vector<double,2> unit_vec_between_nodes(2);
138
139 for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
140 iter != neighbouring_node_indices.end();
141 ++iter)
142 {
143 // The method GetAngleBetweenNodes() returns an angle in the range (-pi,pi]
144 alpha = r_mesh.GetAngleBetweenNodes(index, *iter);
145
146 assert(alpha <= M_PI);
147 assert(alpha > -M_PI);
148
149 if (sin(alpha-theta) > DBL_EPSILON)
150 {
151 // Initialise a zero force vector between neighbouring nodes
152 c_vector<double,2> force_between_nodes = zero_vector<double>(2);
153
154 // Iterate over vector of forces present and add up forces between nodes
155 for (std::vector<boost::shared_ptr<AbstractTwoBodyInteractionForce<2> > >::iterator force_iter = mForceCollection.begin();
156 force_iter != mForceCollection.end();
157 ++force_iter)
158 {
159 force_between_nodes += (*force_iter)->CalculateForceBetweenNodes(index, *iter, mrCellPopulation);
160 }
161
162 unit_vec_between_nodes[0] = cos(alpha);
163 unit_vec_between_nodes[1] = sin(alpha);
164
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);
168 }
169 }
170
171 std::vector<double> ret(2);
172 ret[0] = tangential_force;
173 ret[1] = normal_force;
174
175 return ret;
176}
177
178std::vector<double> DiscreteSystemForceCalculator::GetSamplingAngles(unsigned index)
179{
181
182 std::set<unsigned> neighbouring_node_indices = mrCellPopulation.GetNeighbouringNodeIndices(index);
183
184 std::vector<double> sampling_angles(4*neighbouring_node_indices.size());
185
186 unsigned i=0;
187
188 for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
189 iter != neighbouring_node_indices.end();
190 ++iter)
191 {
192 // The method GetAngleBetweenNodes() returns an angle in the range (-pi,pi]
193 double alpha = r_mesh.GetAngleBetweenNodes(index, *iter);
194
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;
199
200 // Calculate sampling angles in the range (-pi,pi]
201 if (alpha_minus_epsilon <= -M_PI)
202 {
203 alpha_minus_epsilon += 2*M_PI;
204 }
205 sampling_angles[i] = alpha_minus_epsilon;
206
207 assert(sampling_angles[i] <= M_PI);
208 assert(sampling_angles[i] > -M_PI);
209 i++;
210
211 if (alpha_plus_epsilon > M_PI)
212 {
213 alpha_plus_epsilon -= 2*M_PI;
214 }
215 sampling_angles[i] = alpha_plus_epsilon;
216
217 assert(sampling_angles[i] <= M_PI);
218 assert(sampling_angles[i] > -M_PI);
219 i++;
220
221 if (alpha_plus_pi_minus_epsilon > M_PI)
222 {
223 alpha_plus_pi_minus_epsilon -= 2*M_PI;
224 }
225 sampling_angles[i] = alpha_plus_pi_minus_epsilon;
226
227 assert(sampling_angles[i] <= M_PI);
228 assert(sampling_angles[i] > -M_PI);
229 i++;
230
231 if (alpha_plus_pi_plus_epsilon > M_PI)
232 {
233 alpha_plus_pi_plus_epsilon -= 2*M_PI;
234 }
235 sampling_angles[i] = alpha_plus_pi_plus_epsilon;
236
237 assert(sampling_angles[i] <= M_PI);
238 assert(sampling_angles[i] > -M_PI);
239 i++;
240 }
241
242 sort(sampling_angles.begin(), sampling_angles.end());
243 return sampling_angles;
244}
245
246double DiscreteSystemForceCalculator::GetLocalExtremum(unsigned index, double angle1, double angle2)
247{
248 // We always pass in angle1 and angle2 such that angle1<angle2,
249 // but note that angle1 may be <M_PI
250 assert(angle1 < angle2);
251
252 double tolerance = 1e-5;
253 double previous_angle;
254 double current_error;
255 double current_angle = angle1;
256
257 current_error = angle2 - angle1;
258 std::vector<double> current_ft_and_fn(2);
259
260 while (current_error > tolerance)
261 {
262 previous_angle = current_angle;
263 current_ft_and_fn = CalculateFtAndFn(index, current_angle);
264 current_angle -= current_ft_and_fn[0]/current_ft_and_fn[1];
265 current_error = fabs(current_angle - previous_angle);
266 }
267
268 assert(current_angle>angle1 && current_angle<angle2);
269 assert(current_error < tolerance);
270
271 return current_angle;
272}
273
274std::vector<double> DiscreteSystemForceCalculator::GetExtremalAngles(unsigned index, std::vector<double> samplingAngles)
275{
276 std::vector<double> extremal_angles;
277 std::vector<double> ft_and_fn(2);
278 std::vector<double> tangential_force(samplingAngles.size());
279
280 for (unsigned i=0; i<samplingAngles.size(); i++)
281 {
282 ft_and_fn = CalculateFtAndFn(index, samplingAngles[i]);
283 tangential_force[i] = ft_and_fn[0];
284 }
285
286 unsigned n = samplingAngles.size()-1;
287
288 for (unsigned i=0; i<n; i++)
289 {
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))
292 {
293 double next_extremal_angle;
294
295 // If we are in the interval that crosses the branch line at pi,
296 // then subtract 2*pi from the positive angle
297 if (i==n-1)
298 {
299 samplingAngles[i%n] -= 2*M_PI;
300 }
301
302 if (samplingAngles[(i+1)%n] - samplingAngles[i%n] < 2*mEpsilon + 1e-6 )
303 {
304 // If we find a jump through zero, then the local extremum is
305 // simply at the mid-point of the interval
306 next_extremal_angle = 0.5*(samplingAngles[(i+1)%n] + samplingAngles[i%n]);
307 }
308 else
309 {
310 // Otherwise we need to find it using Newton's method
311 next_extremal_angle = GetLocalExtremum(index, samplingAngles[i%n], samplingAngles[(i+1)%n]);
312 }
313
314 if (next_extremal_angle <= -M_PI)
315 {
316 next_extremal_angle += 2*M_PI;
317 }
318 assert(next_extremal_angle>-M_PI && next_extremal_angle<=M_PI);
319 extremal_angles.push_back(next_extremal_angle);
320 }
321 }
322
323 return extremal_angles;
324}
325
327{
328 mEpsilon = epsilon;
329}
virtual unsigned GetNumNodes() const
Node< SPACE_DIM > * GetNode(unsigned index) const
double GetLocalExtremum(unsigned index, double angle1, double angle2)
std::vector< double > GetExtremalAngles(unsigned index, std::vector< double > samplingAngles)
std::vector< boost::shared_ptr< AbstractTwoBodyInteractionForce< 2 > > > mForceCollection
void WriteResultsToFile(std::string simulationOutputDirectory)
std::vector< double > GetSamplingAngles(unsigned index)
DiscreteSystemForceCalculator(MeshBasedCellPopulation< 2 > &rCellPopulation, std::vector< boost::shared_ptr< AbstractTwoBodyInteractionForce< 2 > > > forceCollection)
std::vector< double > CalculateFtAndFn(unsigned index, double theta)
MeshBasedCellPopulation< 2 > & mrCellPopulation
std::vector< std::vector< double > > CalculateExtremalNormalForces()
std::set< unsigned > GetNeighbouringNodeIndices(unsigned index)
MutableMesh< ELEMENT_DIM, SPACE_DIM > & rGetMesh()
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
double GetTime() const
static SimulationTime * Instance()
double GetAngleBetweenNodes(unsigned indexA, unsigned indexB)