cellular_raza_building_blocks/cell_building_blocks/
bacterial_rods.rs

1use crate::{CartesianCuboid, CartesianSubDomain};
2use cellular_raza_concepts::*;
3
4#[cfg(feature = "approx")]
5use approx::AbsDiffEq;
6use num::FromPrimitive;
7use serde::{Deserialize, Serialize};
8
9use nalgebra::{Const, Dyn, Matrix, VecStorage};
10
11/// A mechanical model for Bacterial Rods
12///
13/// See the [Bacterial Rods](https://cellular-raza.com/showcase/bacterial-rods) example for a
14/// detailed example.
15///
16/// # Parameters & Variables
17///
18/// | Symbol | Struct Field | Description |
19/// | --- | --- | --- |
20/// | $\gamma$ | `spring_tension` | Tension of the springs connecting the vertices. |
21/// | $D$ | `diffusion_constant` | Diffusion constant corresponding to brownian motion. |
22/// | $\lambda$ | `damping` | Damping constant. |
23/// | $l$ | `spring_length` | Length of an individual segment between two vertices. |
24/// | $\eta$ | `rigidity` | Rigidity with respect to bending the rod. |
25///
26/// # Equations
27///
28/// The vertices which are being modeled are stored in the `pos` struct field and their
29/// corresponding velocities in the `vel` field.
30///
31/// \\begin{equation}
32///     \vec{v}_i= \text{\texttt{rod\\_mechanics\.pos\.row(i)}}
33/// \\end{equation}
34///
35/// We define the edge $\vec{c}\_i:=\vec{v}\_i-\vec{v}\_{i-1}$.
36/// The first force acts between the vertices $v\_i$ of the model and aims to maintain an equal
37/// distance between all vertices via
38///
39/// \\begin{equation}
40///     \vec{F}\_{i,\text{springs}} = -\gamma\left(1-\frac{l}{||\vec{c}\_i||}\right)\vec{c}\_i
41///         +\gamma\left(1-\frac{l}{||\vec{c}\_{i+1}||}\right)\vec{c}\_{i+1}.
42/// \\end{equation}
43///
44/// We assume the properties of a simple elastic rod.
45/// With the angle $\alpha_i$ between adjacent edges $\vec{c}\_{i-1},\vec{c}\_i$ we can formulate
46/// the bending force which is proportional to the curvature $\kappa\_i$ at vertex $i$
47///
48/// \\begin{equation}
49///     \kappa\_i = 2\tan\left(\frac{\alpha\_i}{2}\right).
50/// \\end{equation}
51///
52/// The resulting force acts along the angle bisector which can be calculated from the edge
53/// vectors.
54/// The forces acting on vertices $\vec{v}\_i,\vec{v}\_{i-1},\vec{v}\_{i+1}$ are given by
55///
56/// \\begin{align}
57///     \vec{F}\_{i,\text{curvature}} &= \eta\kappa_i
58///         \frac{\vec{c}\_i - \vec{c}\_{i+1}}{|\vec{c}\_i-\vec{c}\_{i+1}|}\\\\
59///     \vec{F}\_{i-1,\text{curvature}} &= -\frac{1}{2}\vec{F}\_{i,\text{curvature}}\\\\
60///     \vec{F}\_{i+1,\text{curvature}} &= -\frac{1}{2}\vec{F}\_{i,\text{curvature}}
61/// \\end{align}
62///
63/// where $\eta\_i$ is the angle curvature at vertex $\vec{v}\_i$.
64/// The total force $\vec{F}_{i,\text{total}}$ at vertex $i$ consists of multiple contributions.
65///
66/// \\begin{equation}
67///     \vec{F}\_{i,\text{total}} = \vec{F}\_{i,\text{springs}} + \vec{F}\_{i,\text{curvature}}
68///         + \vec{F}\_{i,\text{external}}
69/// \\end{equation}
70///
71/// # References
72///
73#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
74#[serde(bound = "
75F: 'static
76    + PartialEq
77    + Clone
78    + core::fmt::Debug
79    + Serialize
80    + for<'a> Deserialize<'a>,
81")]
82pub struct RodMechanics<F, const D: usize> {
83    /// The current position
84    pub pos: Matrix<
85        F,
86        nalgebra::Dyn,
87        nalgebra::Const<D>,
88        nalgebra::VecStorage<F, nalgebra::Dyn, nalgebra::Const<D>>,
89    >,
90    /// The current velocity
91    pub vel: Matrix<
92        F,
93        nalgebra::Dyn,
94        nalgebra::Const<D>,
95        nalgebra::VecStorage<F, nalgebra::Dyn, nalgebra::Const<D>>,
96    >,
97    /// Controls magnitude of stochastic motion
98    pub diffusion_constant: F,
99    /// Spring tension between individual vertices
100    pub spring_tension: F,
101    /// Stiffness at each joint connecting two edges
102    pub rigidity: F,
103    /// Target spring length
104    pub spring_length: F,
105    /// Daming constant
106    pub damping: F,
107}
108
109#[cfg(feature = "approx")]
110impl<F, const D: usize> AbsDiffEq for RodMechanics<F, D>
111where
112    F: AbsDiffEq + nalgebra::Scalar,
113    F::Epsilon: Clone,
114{
115    type Epsilon = F::Epsilon;
116
117    fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool {
118        let RodMechanics {
119            pos,
120            vel,
121            diffusion_constant,
122            spring_tension,
123            rigidity,
124            spring_length,
125            damping,
126        } = &self;
127        pos.iter()
128            .zip(other.pos.iter())
129            .all(|(x, y)| x.abs_diff_eq(y, epsilon.clone()))
130            && vel
131                .iter()
132                .zip(other.vel.iter())
133                .all(|(x, y)| x.abs_diff_eq(y, epsilon.clone()))
134            && diffusion_constant.abs_diff_eq(&other.diffusion_constant, epsilon.clone())
135            && spring_tension.abs_diff_eq(&other.spring_tension, epsilon.clone())
136            && rigidity.abs_diff_eq(&other.rigidity, epsilon.clone())
137            && spring_length.abs_diff_eq(&other.spring_length, epsilon.clone())
138            && damping.abs_diff_eq(&other.damping, epsilon)
139    }
140
141    fn default_epsilon() -> Self::Epsilon {
142        F::default_epsilon()
143    }
144}
145
146impl<F, const D: usize>
147    Mechanics<
148        Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
149        Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
150        Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
151        F,
152    > for RodMechanics<F, D>
153where
154    F: nalgebra::RealField + Clone + num::Float,
155    rand_distr::StandardNormal: rand_distr::Distribution<F>,
156{
157    fn calculate_increment(
158        &self,
159        force: Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
160    ) -> Result<
161        (
162            Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
163            Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
164        ),
165        CalcError,
166    > {
167        use core::ops::AddAssign;
168        let one_half = F::one() / (F::one() + F::one());
169        let two = F::one() + F::one();
170
171        // Calculate internal force between the two points of the Agent
172        let mut total_force = force;
173
174        // Calculate force exerted by spring action between individual vertices
175        let n_rows = self.pos.nrows();
176        let dist_internal = self.pos.rows(0, n_rows - 1) - self.pos.rows(1, n_rows - 1);
177        dist_internal.row_iter().enumerate().for_each(|(i, dist)| {
178            if !dist.norm().is_zero() {
179                let dir = dist.normalize();
180                let force_internal =
181                    -dir * self.spring_tension * (dist.norm() - self.spring_length);
182                total_force.row_mut(i).add_assign(force_internal * one_half);
183                total_force
184                    .row_mut(i + 1)
185                    .add_assign(-force_internal * one_half);
186            }
187        });
188
189        // Calculate force exerted by curvature-contributions
190        use itertools::Itertools;
191        dist_internal
192            .row_iter()
193            .tuple_windows::<(_, _)>()
194            .enumerate()
195            .for_each(|(i, (conn1, conn2))| {
196                let angle = conn1.angle(&conn2);
197                let force_d = conn1 - conn2;
198                let force_dir = if !force_d.norm().is_zero() {
199                    force_d.normalize()
200                } else {
201                    force_d
202                };
203                let force =
204                    force_dir * two * self.rigidity * <F as num::Float>::tan(one_half * angle);
205                total_force.row_mut(i).add_assign(-force * one_half);
206                total_force.row_mut(i + 1).add_assign(force);
207                total_force.row_mut(i + 2).add_assign(-force * one_half);
208            });
209
210        // Calculate damping force
211        total_force -= &self.vel * self.damping;
212        Ok((self.vel.clone(), total_force))
213    }
214
215    fn get_random_contribution(
216        &self,
217        rng: &mut rand_chacha::ChaCha8Rng,
218        dt: F,
219    ) -> Result<
220        (
221            Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
222            Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
223        ),
224        RngError,
225    > {
226        let distr = match rand_distr::Normal::new(F::zero(), <F as num::Float>::sqrt(dt)) {
227            Ok(e) => Ok(e),
228            Err(e) => Err(cellular_raza_concepts::RngError(format!("{e}"))),
229        }?;
230        let dpos = nalgebra::Matrix::<F, Dyn, Const<D>, _>::from_distribution(
231            self.pos.nrows(),
232            &distr,
233            rng,
234        ) * <F as num::Float>::powi(F::one() + F::one(), -2)
235            * self.diffusion_constant
236            / dt;
237        let dvel = nalgebra::Matrix::<F, Dyn, Const<D>, _>::zeros(self.pos.nrows());
238
239        Ok((dpos, dvel))
240    }
241}
242
243impl<F: Clone, const D: usize> Position<Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>>
244    for RodMechanics<F, D>
245{
246    fn pos(&self) -> Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>> {
247        self.pos.clone()
248    }
249
250    fn set_pos(&mut self, position: &Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>) {
251        self.pos = position.clone();
252    }
253}
254
255impl<F: Clone, const D: usize> Velocity<Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>>
256    for RodMechanics<F, D>
257{
258    fn velocity(&self) -> Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>> {
259        self.vel.clone()
260    }
261
262    fn set_velocity(&mut self, velocity: &Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>) {
263        self.vel = velocity.clone();
264    }
265}
266
267/// Automatically derives a [Interaction] suitable for rods from a point-wise interaction.
268#[derive(Clone, Debug, Deserialize, PartialEq, Serialize)]
269#[cfg_attr(feature = "approx", derive(AbsDiffEq))]
270pub struct RodInteraction<I>(pub I);
271
272impl<I, Inf> InteractionInformation<Inf> for RodInteraction<I>
273where
274    I: InteractionInformation<Inf>,
275{
276    fn get_interaction_information(&self) -> Inf {
277        self.0.get_interaction_information()
278    }
279}
280
281impl<I, F, Inf, const D: usize>
282    Interaction<
283        Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
284        Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
285        Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
286        Inf,
287    > for RodInteraction<I>
288where
289    I: Interaction<nalgebra::SVector<F, D>, nalgebra::SVector<F, D>, nalgebra::SVector<F, D>, Inf>,
290    F: 'static + nalgebra::RealField + Copy + core::fmt::Debug + num::Zero,
291{
292    fn calculate_force_between(
293        &self,
294        own_pos: &Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
295        own_vel: &Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
296        ext_pos: &Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
297        ext_vel: &Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
298        ext_inf: &Inf,
299    ) -> Result<
300        (
301            Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
302            Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
303        ),
304        CalcError,
305    > {
306        use core::ops::AddAssign;
307        use itertools::Itertools;
308        let mut force_own = nalgebra::Matrix::<F, Dyn, Const<D>, _>::zeros(own_vel.nrows());
309        let mut force_ext = nalgebra::Matrix::<F, Dyn, Const<D>, _>::zeros(own_vel.nrows());
310        for (i, p1) in own_pos.row_iter().enumerate() {
311            for (j, (p2_n0, p2_n1)) in ext_pos.row_iter().tuple_windows::<(_, _)>().enumerate() {
312                // Calculate the closest point of the external position
313                let (_, nearest_point, rel_length) = crate::nearest_point_from_point_to_line(
314                    &p1.transpose(),
315                    &(p2_n0.transpose(), p2_n1.transpose()),
316                );
317
318                let (f_own, f_ext) = self.0.calculate_force_between(
319                    &p1.transpose(),
320                    &own_vel.row(i).transpose(),
321                    &nearest_point,
322                    &ext_vel.row(j).transpose(),
323                    ext_inf,
324                )?;
325
326                force_own.row_mut(i).add_assign(f_own.transpose());
327                force_ext
328                    .row_mut(j)
329                    .add_assign(f_ext.transpose() * (F::one() - rel_length));
330                force_ext
331                    .row_mut((j + 1) % own_pos.nrows())
332                    .add_assign(f_ext.transpose() * rel_length);
333            }
334        }
335        Ok((force_own, force_ext))
336    }
337
338    fn is_neighbor(
339        &self,
340        own_pos: &Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
341        ext_pos: &Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
342        ext_inf: &Inf,
343    ) -> Result<bool, CalcError> {
344        for p in own_pos.row_iter() {
345            for q in ext_pos.row_iter() {
346                if self
347                    .0
348                    .is_neighbor(&p.transpose(), &q.transpose(), ext_inf)?
349                {
350                    return Ok(true);
351                }
352            }
353        }
354        Ok(false)
355    }
356
357    fn react_to_neighbors(&mut self, neighbors: usize) -> Result<(), CalcError> {
358        self.0.react_to_neighbors(neighbors)
359    }
360}
361
362/// Cells are represented by rods
363#[derive(Domain)]
364pub struct CartesianCuboidRods<F, const D: usize> {
365    /// The base-cuboid which is being repurposed
366    #[DomainRngSeed]
367    pub domain: CartesianCuboid<F, D>,
368    /// Gel pressure which is only relevant for 3D simulations and always acts with constant
369    /// force downwards (negative z-direction).
370    pub gel_pressure: F,
371    /// Computes friction at all surfaces of the box
372    pub surface_friction: F,
373    /// The distance for which the friction will be applied
374    pub surface_friction_distance: F,
375}
376
377impl<C, F, const D: usize> SortCells<C> for CartesianCuboidRods<F, D>
378where
379    C: Position<Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>>,
380    F: 'static
381        + nalgebra::Field
382        + Clone
383        + core::fmt::Debug
384        + num::FromPrimitive
385        + num::ToPrimitive
386        + num::Float
387        + Copy,
388{
389    type VoxelIndex = [usize; D];
390
391    fn get_voxel_index_of(&self, cell: &C) -> Result<Self::VoxelIndex, BoundaryError> {
392        let pos = cell.pos().row_sum().transpose() / F::from_usize(cell.pos().nrows()).unwrap();
393        let index = self.domain.get_voxel_index_of_raw(&pos)?;
394        Ok(index)
395    }
396}
397
398impl<F, const D: usize> DomainCreateSubDomains<CartesianSubDomainRods<F, D>>
399    for CartesianCuboidRods<F, D>
400where
401    F: 'static + num::Float + core::fmt::Debug + num::FromPrimitive,
402{
403    type SubDomainIndex = usize;
404    type VoxelIndex = [usize; D];
405
406    fn create_subdomains(
407        &self,
408        n_subdomains: std::num::NonZeroUsize,
409    ) -> Result<
410        impl IntoIterator<
411            Item = (
412                Self::SubDomainIndex,
413                CartesianSubDomainRods<F, D>,
414                Vec<Self::VoxelIndex>,
415            ),
416        >,
417        DecomposeError,
418    > {
419        let subdomains = self.domain.create_subdomains(n_subdomains)?;
420        Ok(subdomains
421            .into_iter()
422            .map(move |(subdomain_index, subdomain, voxels)| {
423                (
424                    subdomain_index,
425                    CartesianSubDomainRods::<F, D> {
426                        subdomain,
427                        gel_pressure: self.gel_pressure,
428                        surface_friction: self.surface_friction,
429                        surface_friction_distance: self.surface_friction,
430                    },
431                    voxels,
432                )
433            }))
434    }
435}
436
437/// The corresponding SubDomain of the [CartesianCuboidRods] domain.
438#[derive(Clone, SubDomain, Serialize, Deserialize)]
439#[serde(bound = "
440F: 'static
441    + PartialEq
442    + Clone
443    + core::fmt::Debug
444    + Serialize
445    + for<'a> Deserialize<'a>,
446[usize; D]: Serialize + for<'a> Deserialize<'a>,
447")]
448pub struct CartesianSubDomainRods<F, const D: usize> {
449    /// Base subdomain as created by the [CartesianCuboid] domain.
450    #[Base]
451    pub subdomain: CartesianSubDomain<F, D>,
452    /// See [CartesianCuboidRods]
453    pub gel_pressure: F,
454    /// Computes friction at all surfaces of the box
455    pub surface_friction: F,
456    /// The distance for which the friction will be applied
457    pub surface_friction_distance: F,
458}
459
460impl<F>
461    SubDomainForce<
462        Matrix<F, Dyn, Const<3>, VecStorage<F, Dyn, Const<3>>>,
463        Matrix<F, Dyn, Const<3>, VecStorage<F, Dyn, Const<3>>>,
464        Matrix<F, Dyn, Const<3>, VecStorage<F, Dyn, Const<3>>>,
465        F,
466    > for CartesianSubDomainRods<F, 3>
467where
468    F: nalgebra::RealField + num::Float,
469{
470    fn calculate_custom_force(
471        &self,
472        pos: &Matrix<F, Dyn, Const<3>, VecStorage<F, Dyn, Const<3>>>,
473        vel: &Matrix<F, Dyn, Const<3>, VecStorage<F, Dyn, Const<3>>>,
474        _: &F,
475    ) -> Result<
476        Matrix<F, Dyn, Const<3>, VecStorage<F, Dyn, Const<3>>>,
477        cellular_raza_concepts::CalcError,
478    > {
479        use core::ops::AddAssign;
480        let mut force = nalgebra::MatrixXx3::from_fn(pos.nrows(), |_, m| {
481            if m == 2 {
482                -self.gel_pressure
483            } else {
484                F::zero()
485            }
486        });
487        for (i, (p, v)) in pos.row_iter().zip(vel.row_iter()).enumerate() {
488            let d1 = (p.transpose() - self.subdomain.domain_min)
489                .map(|x| <F as num::Float>::abs(x) <= self.surface_friction_distance);
490            let d2 = (p.transpose() - self.subdomain.domain_max)
491                .map(|x| <F as num::Float>::abs(x) <= self.surface_friction_distance);
492            let q = v.norm();
493            if q != F::zero() && d1.zip_map(&d2, |x, y| x || y).into_iter().any(|x| *x) {
494                let dir = v / q;
495                force
496                    .row_mut(i)
497                    .add_assign(-dir * self.gel_pressure * self.surface_friction);
498            }
499        }
500        Ok(force)
501    }
502}
503
504impl<F, const D: usize>
505    SubDomainMechanics<
506        Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
507        Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
508    > for CartesianSubDomainRods<F, D>
509where
510    F: nalgebra::RealField + num::Float,
511{
512    fn apply_boundary(
513        &self,
514        pos: &mut Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
515        vel: &mut Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>,
516    ) -> Result<(), BoundaryError> {
517        // TODO refactor this with matrix multiplication!!!
518        // This will probably be much more efficient and less error-prone!
519
520        // For each position in the springs Agent<D1, D2>
521        let two = F::one() + F::one();
522        pos.row_iter_mut()
523            .zip(vel.row_iter_mut())
524            .for_each(|(mut p, mut v)| {
525                // For each dimension in the space
526                for i in 0..p.ncols() {
527                    // Check if the particle is below lower edge
528                    if p[i] < self.subdomain.get_domain_min()[i] {
529                        p[i] = self.subdomain.get_domain_min()[i] * two - p[i];
530                        v[i] = <F as num::Float>::abs(v[i]);
531                    }
532
533                    // Check if the particle is over the edge
534                    if p[i] > self.subdomain.get_domain_max()[i] {
535                        p[i] = self.subdomain.get_domain_max()[i] * two - p[i];
536                        v[i] = -<F as num::Float>::abs(v[i]);
537                    }
538                }
539            });
540
541        // If new pos is still out of boundary return error
542        for j in 0..pos.nrows() {
543            let p = pos.row(j);
544            for i in 0..pos.ncols() {
545                if p[i] < self.subdomain.get_domain_min()[i]
546                    || p[i] > self.subdomain.get_domain_max()[i]
547                {
548                    return Err(BoundaryError(format!(
549                        "Particle is out of domain at pos {:?}",
550                        pos
551                    )));
552                }
553            }
554        }
555        Ok(())
556    }
557}
558
559impl<C, F, const D: usize> SortCells<C> for CartesianSubDomainRods<F, D>
560where
561    C: Position<Matrix<F, Dyn, Const<D>, VecStorage<F, Dyn, Const<D>>>>,
562    F: 'static
563        + nalgebra::Field
564        + Clone
565        + core::fmt::Debug
566        + num::FromPrimitive
567        + num::ToPrimitive
568        + num::Float
569        + Copy,
570{
571    type VoxelIndex = [usize; D];
572
573    fn get_voxel_index_of(&self, cell: &C) -> Result<Self::VoxelIndex, BoundaryError> {
574        let pos = cell.pos().row_sum().transpose() / F::from_usize(cell.pos().nrows()).unwrap();
575        let index = self.subdomain.get_index_of(pos)?;
576        Ok(index)
577    }
578}
579
580impl<F, const D: usize> RodMechanics<F, D> {
581    /// Divides a [RodMechanics] struct into two thus separating their positions
582    ///
583    /// ```
584    /// # use cellular_raza_building_blocks::*;
585    /// use nalgebra::MatrixXx2;
586    /// let n_vertices = 7;
587    /// let mut pos = MatrixXx2::zeros(n_vertices);
588    /// pos
589    ///     .row_iter_mut()
590    ///     .enumerate()
591    ///     .for_each(|(n_row, mut r)| r[0] += n_row as f32 * 0.5);
592    /// let mut m1 = RodMechanics {
593    ///     pos,
594    ///     vel: MatrixXx2::zeros(n_vertices),
595    ///     diffusion_constant: 0.0,
596    ///     spring_tension: 0.1,
597    ///     rigidity: 0.05,
598    ///     spring_length: 0.5,
599    ///     damping: 0.0,
600    /// };
601    /// let radius = 0.25;
602    /// let m2 = m1.divide(radius)?;
603    ///
604    /// let last_pos_m1 = m1.pos.row(6);
605    /// let first_pos_m2 = m2.pos.row(0);
606    /// assert!(((last_pos_m1 - first_pos_m2).norm() - 2.0 * radius).abs() < 1e-3);
607    /// # Result::<(), cellular_raza_concepts::DivisionError>::Ok(())
608    /// ```
609    pub fn divide(
610        &mut self,
611        radius: F,
612    ) -> Result<RodMechanics<F, D>, cellular_raza_concepts::DivisionError>
613    where
614        F: num::Float + nalgebra::RealField + FromPrimitive + std::iter::Sum,
615    {
616        use itertools::Itertools;
617        let pos = self.pos();
618        let c1 = self;
619        let mut c2 = c1.clone();
620
621        let n_rows = c1.pos.nrows();
622        // Calculate the fraction of how much we need to scale down the individual spring lengths
623        // in order for the distances to still work.
624        let two = F::one() + F::one();
625        let one_half = F::one() / two;
626        let div_factor = one_half - radius / (F::from_usize(n_rows).unwrap() * c1.spring_length);
627
628        // Shrink spring length
629        let new_spring_length = div_factor * c1.spring_length;
630        c1.spring_length = new_spring_length;
631        c2.spring_length = new_spring_length;
632
633        // Set starting point
634        c1.pos.set_row(0, &pos.row(0));
635        c2.pos
636            .set_row(c2.pos.nrows() - 1, &pos.row(c2.pos.nrows() - 1));
637
638        let segments: Vec<_> = pos
639            .row_iter()
640            .tuple_windows::<(_, _)>()
641            .map(|(x, y)| (x - y).norm())
642            .collect();
643        let segment_length = (segments.iter().copied().sum::<F>() - two * radius)
644            / F::from_usize(c2.pos.nrows() - 1).unwrap()
645            / two;
646
647        for n_vertex in 0..c2.pos.nrows() {
648            // Get smallest index k such that the beginning of the new segment is "farther" than the
649            // original vertex at this index k.
650            let k = (0..segments.len())
651                .filter(|n| {
652                    segments.iter().copied().take(*n).sum::<F>()
653                        <= F::from_usize(n_vertex).unwrap() * segment_length
654                })
655                .max()
656                .ok_or(DivisionError(
657                    "Could not calculate new segment lengths".to_string(),
658                ))?;
659            let q = (F::from_usize(n_vertex).unwrap() * segment_length
660                - segments.iter().copied().take(k).sum::<F>())
661                / segments[k];
662            c1.pos.set_row(
663                n_vertex,
664                &(pos.row(k) * (F::one() - q) + pos.row(k + 1) * q),
665            );
666
667            let m = (0..c2.pos.nrows())
668                .filter(|n| {
669                    segments.iter().rev().copied().take(*n).sum::<F>()
670                        <= F::from_usize(n_vertex).unwrap() * segment_length
671                })
672                .max()
673                .ok_or(DivisionError(
674                    "Could not calculate new segment lengths".to_string(),
675                ))?;
676            let p = (F::from_usize(n_vertex).unwrap() * segment_length
677                - segments.iter().rev().copied().take(m).sum::<F>())
678                / segments[c2.pos.nrows() - m - 2];
679            c2.pos.set_row(
680                c2.pos.nrows() - n_vertex - 1,
681                &(pos.row(c2.pos.nrows() - m - 1) * (F::one() - p)
682                    + pos.row(c2.pos.nrows() - m - 2) * p),
683            );
684        }
685        Ok(c2)
686    }
687}
cellular_raza_building_blocks/cell_building_blocks/bacterial_rods.rs

cellular_raza_building_blocks/cell_building_blocks/
bacterial_rods.rs
