BioDynaMo  v1.03.61-fe34ced0
interaction_force.cc
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 CERN & University of Surrey for the benefit of the
4 // BioDynaMo collaboration. All Rights Reserved.
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 //
9 // See the LICENSE file distributed with this work for details.
10 // See the NOTICE file distributed with this work for additional information
11 // regarding copyright ownership.
12 //
13 // -----------------------------------------------------------------------------
14 
15 #include "core/interaction_force.h"
16 
17 #include <algorithm>
18 #include <cmath>
19 
20 #include "core/agent/agent.h"
21 #include "core/shape.h"
22 #include "core/simulation.h"
23 #include "core/util/log.h"
24 #include "core/util/math.h"
25 #include "core/util/random.h"
26 #include "core/util/type.h"
28 
29 namespace bdm {
30 
31 using neuroscience::NeuriteElement;
32 
33 Real4 InteractionForce::Calculate(const Agent* lhs, const Agent* rhs) const {
34  if (lhs->GetShape() == Shape::kSphere && rhs->GetShape() == Shape::kSphere) {
35  Real3 result;
36  ForceBetweenSpheres(lhs, rhs, &result);
37  return {result[0], result[1], result[2], 0};
38  } else if (lhs->GetShape() == Shape::kSphere &&
39  rhs->GetShape() == Shape::kCylinder) {
40  Real3 result;
41  ForceOnASphereFromACylinder(lhs, rhs, &result);
42  return {result[0], result[1], result[2], 0};
43  } else if (lhs->GetShape() == Shape::kCylinder &&
44  rhs->GetShape() == Shape::kSphere) {
45  Real4 result;
46  ForceOnACylinderFromASphere(lhs, rhs, &result);
47  return result;
48  } else if (lhs->GetShape() == Shape::kCylinder &&
49  rhs->GetShape() == Shape::kCylinder) {
50  Real4 result;
51  ForceBetweenCylinders(lhs, rhs, &result);
52  return result;
53  } else {
54  Log::Fatal("InteractionForce",
55  "InteractionForce only supports sphere or cylinder shapes");
56  return {0, 0, 0, 0};
57  }
58 }
59 
61  const Agent* sphere_rhs,
62  Real3* result) const {
63  const Real3& ref_mass_location = sphere_lhs->GetPosition();
64  real_t ref_diameter = sphere_lhs->GetDiameter();
65  real_t ref_iof_coefficient = 0.15;
66  const Real3& nb_mass_location = sphere_rhs->GetPosition();
67  real_t nb_diameter = sphere_rhs->GetDiameter();
68  real_t nb_iof_coefficient = 0.15;
69 
70  auto c1 = ref_mass_location;
71  real_t r1 = 0.5 * ref_diameter;
72  auto c2 = nb_mass_location;
73  real_t r2 = 0.5 * nb_diameter;
74  // We take virtual bigger radii to have a distant interaction, to get a
75  // desired density.
76  real_t additional_radius =
77  10.0 * std::min(ref_iof_coefficient, nb_iof_coefficient);
78  r1 += additional_radius;
79  r2 += additional_radius;
80  // the 3 components of the vector c2 -> c1
81  real_t comp1 = c1[0] - c2[0];
82  real_t comp2 = c1[1] - c2[1];
83  real_t comp3 = c1[2] - c2[2];
84  real_t center_distance =
85  std::sqrt(comp1 * comp1 + comp2 * comp2 + comp3 * comp3);
86  // the overlap distance (how much one penetrates in the other)
87  real_t delta = r1 + r2 - center_distance;
88  // if no overlap : no force
89  if (delta < 0) {
90  *result = {0.0, 0.0, 0.0};
91  return;
92  }
93  // to avoid a division by 0 if the centers are (almost) at the same
94  // location
95  if (center_distance < 0.00000001) {
96  auto* random = Simulation::GetActive()->GetRandom();
97  auto force2on1 = random->template UniformArray<3>(-3.0, 3.0);
98  *result = force2on1;
99  return;
100  }
101  // the force itself
102  real_t r = (r1 * r2) / (r1 + r2);
103  real_t gamma = 1; // attraction coeff
104  real_t k = 2; // repulsion coeff
105  real_t f = k * delta - gamma * std::sqrt(r * delta);
106 
107  real_t force_module = f / center_distance;
108  Real3 force2on1(
109  {force_module * comp1, force_module * comp2, force_module * comp3});
110  *result = force2on1;
111 }
112 
114  const Agent* sphere,
115  Real4* result) const {
116  auto* ne = bdm_static_cast<const NeuriteElement*>(cylinder);
117  auto proximal_end = ne->ProximalEnd();
118  auto distal_end = ne->DistalEnd();
119  auto axis = ne->GetSpringAxis();
120  // TODO(neurites) use cylinder.GetActualLength() ??
121  real_t actual_length = axis.Norm();
122  real_t d = ne->GetDiameter();
123  auto c = sphere->GetPosition();
124  real_t r = 0.5 * sphere->GetDiameter();
125 
126  // I. If the cylinder is small with respect to the sphere:
127  // we only consider the interaction between the sphere and the point mass
128  // (i.e. distal point) of the cylinder - that we treat as a sphere.
129  if (actual_length < r) {
130  // move back sphere center by 1 cylinder radius from distal_end
131  // vector_x = rc * (axis[0]/actual_length)
132  // vector_y = rc * (axis[1]/actual_length)
133  // vector_z = rc * (axis[2]/actual_length)
134  real_t rc = 0.5 * d;
135  Real3 dvec = (axis / actual_length) * rc; // displacement vector
136  Real3 npd = distal_end - dvec; // new sphere center
137  *result = ComputeForceOfASphereOnASphere(npd, rc, c, r);
138  return;
139  }
140 
141  // II. If the cylinder is of the same scale or bigger than the sphere,
142  // we look at the interaction between the sphere and the closest point
143  // (to the sphere center) on the cylinder. This interaction is distributed
144  // to
145  // the two ends of the cylinder: the distal (point mass of the segment) and
146  // the proximal (point mass of the mother of the segment).
147 
148  // 1) Finding cc : the closest point to c on the line proximal_end
149  // proximal_end
150  // ("line" and not "segment")
151  // It is the projection of the vector proximal_end->c onto the vector
152  // proximal_end->distal_end
153  // (=axis)
154  auto proximal_end_closest = c - proximal_end;
155 
156  // projection of proximal_end_closest onto axis =
157  // (proximal_end_closest.axis)/norm(axis)^2 * axis
158  // length of the projection = (proximal_end_closest.axis)/norm(axis)
159  real_t proximal_end_closest_axis =
160  proximal_end_closest.EntryWiseProduct(axis).Sum();
161  real_t k = proximal_end_closest_axis / (actual_length * actual_length);
162  // cc = proximal_end + k* axis
163  Real3 cc = proximal_end + (axis * k);
164 
165  // 2) Look if c -and hence cc- is (a) between proximal_end and distal_end,
166  // (b) before proximal_end or
167  // (c) after distal_end
168  real_t proportion_to_proximal_end;
169  if (k <= 1.0 && k >= 0.0) {
170  // a) if cc (the closest point to c on the line pPpD) is between
171  // proximal_end
172  // and distal_end
173  // the force is distributed to the two nodes
174  proportion_to_proximal_end = 1.0 - k;
175  } else if (k < 0) {
176  // b) if the closest point to c on the line pPpD is before
177  // proximal_end
178  // the force is only on the proximal end (the mother point mass)
179  proportion_to_proximal_end = 1.0;
180  cc = proximal_end;
181  } else { // if(k>1)
182  // c) if cc is after distal_end, the force is only on the distal end
183  // (the
184  // segment's point mass).
185  proportion_to_proximal_end = 0.0;
186  cc = distal_end;
187  }
188 
189  // 3) If the smallest distance between the cylinder and the center of the
190  // sphere
191  // is larger than the radius of the two objects , there is no
192  // interaction:
193  real_t penetration = d / 2 + r - Math::GetL2Distance(c, cc);
194  if (penetration <= 0) {
195  *result = Real4{0.0, 0.0, 0.0, 0.0};
196  return;
197  }
198  auto force = ComputeForceOfASphereOnASphere(cc, d * 0.5, c, r);
199  *result = {force[0], force[1], force[2], proportion_to_proximal_end};
200  return;
201 }
202 
204  const Agent* cylinder,
205  Real3* result) const {
206  // it is the opposite of force on a cylinder from sphere:
207  Real4 temp;
208  ForceOnACylinderFromASphere(cylinder, sphere, &temp);
209 
210  *result = {-temp[0], -temp[1], -temp[2]};
211 }
212 
214  const Agent* cylinder2,
215  Real4* result) const {
216  auto* c1 = bdm_static_cast<const NeuriteElement*>(cylinder1);
217  auto* c2 = bdm_static_cast<const NeuriteElement*>(cylinder2);
218  auto a = c1->ProximalEnd();
219  auto b = c1->GetMassLocation();
220  real_t d1 = c1->GetDiameter();
221  auto c = c2->ProximalEnd();
222  auto d = c2->GetMassLocation();
223  real_t d2 = c2->GetDiameter();
224 
225  real_t k = 0.5; // part devoted to the distal node
226 
227  // looking for closest point on them
228  // (based on http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline3d/)
229  real_t p13x = a[0] - c[0];
230  real_t p13y = a[1] - c[1];
231  real_t p13z = a[2] - c[2];
232  real_t p43x = d[0] - c[0];
233  real_t p43y = d[1] - c[1];
234  real_t p43z = d[2] - c[2];
235  real_t p21x = b[0] - a[0];
236  real_t p21y = b[1] - a[1];
237  real_t p21z = b[2] - a[2];
238 
239  real_t d1343 = p13x * p43x + p13y * p43y + p13z * p43z;
240  real_t d4321 = p21x * p43x + p21y * p43y + p21z * p43z;
241  real_t d1321 = p21x * p13x + p21y * p13y + p21z * p13z;
242  real_t d4343 = p43x * p43x + p43y * p43y + p43z * p43z;
243  real_t d2121 = p21x * p21x + p21y * p21y + p21z * p21z;
244 
245  Real3 p1, p2;
246 
247  real_t denom = d2121 * d4343 - d4321 * d4321;
248 
249  // if the two segments are not ABSOLUTLY parallel
250  if (denom > 0.000000000001) {
251  real_t numer = d1343 * d4321 - d1321 * d4343;
252 
253  real_t mua = numer / denom;
254  real_t mub = (d1343 + mua * d4321) / d4343;
255 
256  if (mua < 0) {
257  p1 = a;
258  k = 1;
259  } else if (mua > 1) {
260  p1 = b;
261  k = 0;
262  } else {
263  p1 = Real3{a[0] + mua * p21x, a[1] + mua * p21y, a[2] + mua * p21z};
264  k = 1 - mua;
265  }
266 
267  if (mub < 0) {
268  p2 = c;
269  } else if (mub > 1) {
270  p2 = d;
271  } else {
272  p2 = Real3{c[0] + mub * p43x, c[1] + mub * p43y, c[2] + mub * p43z};
273  }
274 
275  } else {
276  p1 = a + (b - a) * 0.5;
277  p2 = c + (d - c) * 0.5;
278  }
279 
280  // W put a virtual sphere on the two cylinders
281  auto force = ComputeForceOfASphereOnASphere(p1, d1 / 2.0, p2, d2 / 2.0) * 10;
282 
283  *result = {force[0], force[1], force[2], k};
284 }
285 
287  real_t r1,
288  const Real3& c2,
289  real_t r2) const {
290  real_t comp1 = c1[0] - c2[0];
291  real_t comp2 = c1[1] - c2[1];
292  real_t comp3 = c1[2] - c2[2];
293  real_t distance_between_centers =
294  std::sqrt(comp1 * comp1 + comp2 * comp2 + comp3 * comp3);
295  // the overlap distance (how much one penetrates in the other)
296  real_t a = r1 + r2 - distance_between_centers;
297  // if no overlap: no force
298  if (a < 0) {
299  return Real4{0.0, 0.0, 0.0, 0.0};
300  }
301  // to avoid a division by 0 if the centers are (almost) at the same location
302  if (distance_between_centers <
303  0.00000001) { // TODO(neurites) hard coded values
304  auto* random = Simulation::GetActive()->GetRandom();
305  auto force2on1 = random->template UniformArray<3>(-3.0, 3.0);
306  return Real4{force2on1[0], force2on1[1], force2on1[2], 0.0};
307  } else {
308  // the force is prop to the square of the interpentration distance and to
309  // the radii.
310  real_t force_module = a / distance_between_centers;
311  Real4 force2on1({force_module * comp1, force_module * comp2,
312  force_module * comp3, 0.0});
313  return force2on1;
314  }
315 }
316 
317 } // namespace bdm
shape.h
bdm::InteractionForce::Calculate
virtual Real4 Calculate(const Agent *lhs, const Agent *rhs) const
Definition: interaction_force.cc:33
bdm::InteractionForce::ComputeForceOfASphereOnASphere
Real4 ComputeForceOfASphereOnASphere(const Real3 &c1, real_t r1, const Real3 &c2, real_t r2) const
Definition: interaction_force.cc:286
bdm
Definition: agent.cc:39
bdm::real_t
double real_t
Definition: real_t.h:21
bdm::Simulation::GetRandom
Random * GetRandom()
Returns a thread local random number generator.
Definition: simulation.cc:284
bdm::Agent::GetDiameter
virtual real_t GetDiameter() const =0
neurite_element.h
type.h
bdm::Agent::GetShape
virtual Shape GetShape() const =0
bdm::InteractionForce::ForceOnASphereFromACylinder
void ForceOnASphereFromACylinder(const Agent *sphere, const Agent *cylinder, Real3 *result) const
Definition: interaction_force.cc:203
math.h
random.h
bdm::Agent
Contains code required by all agents.
Definition: agent.h:79
bdm::InteractionForce::ForceBetweenCylinders
void ForceBetweenCylinders(const Agent *cylinder1, const Agent *cylinder2, Real4 *result) const
Definition: interaction_force.cc:213
bdm::InteractionForce::ForceBetweenSpheres
void ForceBetweenSpheres(const Agent *sphere_lhs, const Agent *sphere_rhs, Real3 *result) const
Definition: interaction_force.cc:60
bdm::kSphere
@ kSphere
Definition: shape.h:20
agent.h
bdm::Log::Fatal
static void Fatal(const std::string &location, const Args &... parts)
Prints fatal error message.
Definition: log.h:115
log.h
bdm::InteractionForce::ForceOnACylinderFromASphere
void ForceOnACylinderFromASphere(const Agent *cylinder, const Agent *sphere, Real4 *result) const
Definition: interaction_force.cc:113
simulation.h
bdm::Math::GetL2Distance
static real_t GetL2Distance(const Real3 &pos1, const Real3 &pos2)
Definition: math.h:41
bdm::MathArray
Definition: math_array.h:36
bdm::Agent::GetPosition
virtual const Real3 & GetPosition() const =0
bdm::kCylinder
@ kCylinder
Definition: shape.h:20
bdm::Simulation::GetActive
static Simulation * GetActive()
This function returns the currently active Simulation simulation.
Definition: simulation.cc:67
interaction_force.h