Cytosim PI Cytoskeleton Simulator

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- Class List
- Class Index
- Class Hierarchy
- Class Members

List of all members | Public Member Functions | Public Attributes
Meca Class Reference

## Detailed Description

Meca solves the motion of point-objects defined by the dynamic equation:

d vPTS/dt = mobility * mP * ( vBAS + ( mB + mC + mR + mdiffP ) * vPTS )

This equation includes internal forces (bending elasticity for RigidFiber), and all the interactions between the Mecalbe in the system. It is solved in the presence of Brownian motion, and at low Reynolds number ( mass/inertia effects are negligible ).

- Vector vPTS contains all the Mecable coordinates (x, y, z): Fiber, Sphere, Solid and other PointSet. size(vPTS) = DIM * nbPts.
- Vector vBAS is of same size as vPTS, and includes the constant part obtained by linearization of the forces. It includes for instance the positions of Single, calibrated random forces simulating Brownian motion, and also offsets for periodic boundary conditions.
- Matrix mB is the isotropic part obtained after linearization of the forces. It operates similarly and independently on the different dimension X, Y and Z. mB is square of size nbPts, symmetric and sparse.
- Matrix mC is the non-isotropic part obtained after linearization of the forces. mC is square of size DIM*nbPts, symmetric and sparse.

Typically, mB and mC will inherit the stiffness coefficients of the interactions, while vBAS will get forces (stiffness * position). They are set by using interLink(), interLongLink(), interSideLink(), interSlidingLink(), etc.

- mR add the bending elasticity for RigidFiber, or other internal forces. mR is symmetric of size DIM*nbPts, diagonal by blocks, each block corresponding to a Fiber.
- mP applies the projection due to constrained dynamics. For RigidFiber, this maintains the distance between neighboring points (longitudinal incompressibility). mP is symmetric of size DIM*nbPts, diagonal by blocks, each block corresponding to a Fiber. mP is not actually calculated as a matrix: its application on each block is done by Mecable::setSpeedsFromForces()
- mdiffP is a term coming from the derivative of the projection P, and it is added here because it provides better numerical stability in some situations. You can undefine PROJECTION_DIFF in meca.cc to remove mdiffP.

Note: All Links are disabled if the given PointExacts or PointInterpolated have a point in common, because the matrix elements are not calcuated correctly in that case. Generally, such interactions are anyway not desirable. It would correspond for example to a link inside a segment, or between two successive segments on the same Fiber.

## Public Member Functions
real & | base (index_type ix)
| base for force
Vector | pos (index_type ix) const
| position of a point
| Meca ()
| constructor
void | clear ()
| Clear list of Mecable.
void | add (Mecable *ob)
| Add a Mecable to the list of objects to be simulated. More...
unsigned | nbMecables () const
| Number of Mecable.
unsigned | nbPoints () const
| number of points in the system
unsigned | size () const
| Implementation of Solver::LinearOperator.
void | multiply (const real *X, real *Y) const
| calculate Y = M*X, where M is the matrix associated with the system More...
void | precondition (const real *X, real *Y) const
| apply preconditionner: Y = P*X
void | addPureForce (PointExact const &, Vector const &force)
| Add a constant force at a PointExact. More...
void | addPureForce (PointInterpolated const &, Vector const &force)
| Add a constant force at a PointInterpolated. More...
void | addPureTorque (PointInterpolated const &, Torque const &torque)
| Add a torque to a segment.
void | addTorqueClamp (PointInterpolated const &, Vector const &dir, real weight)
| Add a torque to constrain the segment in given direction.
void | interLink (PointExact const &, PointExact const &, real weight)
| Force of stiffness weight. More...
void | interLink (PointInterpolated const &, PointExact const &, real weight)
| Force of stiffness weight. More...
void | interLink (PointInterpolated const &, PointInterpolated const &, real weight)
| Force of stiffness weight. More...
void | interLongLink (PointExact const &, PointExact const &, real len, real weight)
| Force of stiffness weight and resting length len. More...
void | interLongLink (PointInterpolated const &, PointExact const &, real len, real weight)
| Force of stiffness weight and resting length len. More...
void | interLongLink (PointInterpolated const &, PointInterpolated const &, real len, real weight)
| Force of stiffness weight and resting length len. More...
void | interSideLink (PointInterpolated const &, PointExact const &, real len, real weight)
| Force of stiffness weight and resting length len, on the side of first fiber. More...
void | interSideLink (PointInterpolated const &, PointInterpolated const &, real len, real weight)
| Force of stiffness weight and resting length len, on the side of first fiber. More...
void | interSideSideLink (PointInterpolated const &, PointInterpolated const &, real len, real weight)
| Force of stiffness weight and resting length len, on the sides of both fibers. More...
void | interSlidingLink (PointInterpolated const &, PointExact const &, real weight)
| Force of stiffness weight and resting length len, which can slide on first fiber. More...
void | interSlidingLink (PointInterpolated const &, PointInterpolated const &, real weight)
| Force of stiffness weight and resting length len, which can slide on first fiber. More...
void | interSideSlidingLink (PointInterpolated const &, PointExact const &, real len, real weight)
| Force of stiffness weight and resting length len, on the side of first point, which can slide. More...
void | interSideSlidingLink (PointInterpolated const &, PointInterpolated const &, real len, real weight)
| Force of stiffness weight and resting length len, on the side of first point, which can slide. More...
void | interClamp (PointExact const &, const real g[], real weight)
| Force of stiffness weight with fixed position g. More...
void | interClamp (PointInterpolated const &, const real g[], real weight)
| Force of stiffness weight with fixed position g. More...
void | interLongClamp (PointExact const &, Vector const &center, real len, real weight)
| Force of stiffness weight with fixed position g, on the side. More...
void | interSlidingClamp (PointInterpolated const &, Vector const &g, real weight)
| Force of stiffness weight with fixed position g, which can slide. More...
void | interSideClamp (PointInterpolated const &, const real g[], real len, real weight)
| Force of stiffness weight with fixed position g, on the side. More...
void | interPlane (PointExact const &, Vector const &dir, Vector const &g, real weight)
| Force of stiffness weight with a plane defined by g and its normal dir. More...
void | interPlane (PointInterpolated const &, Vector const &dir, Vector const &g, real weight)
| Force of stiffness weight with a plane defined by g and its normal dir. More...
void | interCoulomb (PointExact const &, PointExact const &, real weight)
| Linearized Coulomb repulsive force (experimental)
void | prepare (SimulProp const *)
| Allocate the memory necessary to solve(). This must be called after the last add() More...
void | solve (SimulProp const *, bool precondition)
| Calculate motion of the system. More...
void | calculateForces ()
| calculate Forces on objects and Lagrange multipliers for Fiber, without thermal motion More...
void | printMatrix (std::ostream &) const
| Print the complete matrix (for debugging) More...
void | dumpMatrix (FILE *) const
| Save complete matrix in binary file. More...
void | dumpProjection (FILE *, Mecable const *, bool diff) const
| Save complete matrix in binary file. More...
void | dumpDiagonal (FILE *) const
| Save complete matrix in binary file.
void | dump () const
| Output vectors and matrices in various files (for debugging) More...

## Public Attributes
MatrixSparseSymmetric1 | mB
| isotropic symmetric part of the dynamic, size (nbPts)^2 More...
MatrixSparseSymmetric1 | mC
| non-isotropic symmetric part of the dynamic, size (DIM*nbPts)^2 More...

## Member Function Documentation

void add | ( | Mecable * | o | ) |

Register a new Mecable, set Mecable::matIndex() and update Meca::nbPts and Meca::largestBlock

void addPureForce | ( | PointExact const & | pte,
|  | Vector const & | force
| ) |  |

Display interactions from the Meca Add constant force to pte

void addPureForce | ( | PointInterpolated const & | pti,
|  | Vector const & | force
| ) |  |

Add constant force to pti

void calculateForces | ( |  | ) |

Calculates the force in the objects, that can be accessed by Mecable::netForceP() and calculate the speed of the objects in vRHS, in the abscence of Thermal motion, ie. the motion is purely due to external forces.

This also sets the Lagrange multipliers for the Fiber.

The function will not change the position of the Mecables.

void dump | ( |  | ) | const

This dump the matrix and some vectors in binary files.

Here is some matlab code to read the output:

dim = load('ddim.txt');
mat = fread(fopen('dmat.bin'), [dim, dim], 'double');
rhs = fread(fopen('drhs.bin'), dim, 'double');
sol = fread(fopen('dsol.bin'), dim, 'double');

And to compare the results with matlab's own method, using a scatter plot:

x = bicgstab(mat, rhs, 0.001, 100);
plot(x, sol, '.');

void dumpMatrix | ( | FILE * | file | ) | const

Save the full matrix associated with matVect, for debugging purposes

void dumpProjection | ( | FILE * | file,
|  | Mecable const * | mec,
|  | bool | diff
| ) |  | const

Save the projection matrix

void interClamp | ( | PointExact const & | pta,
|  | const real | g[],
|  | real | weight
| ) |  |

Update Meca to include a link between a point A and a fixed position G The force is linear: force_A = weight * ( G - A ); There is no counter-force in G, since G is immobile.

void interClamp | ( | PointInterpolated const & | pti,
|  | const real | g[],
|  | real | weight
| ) |  |

Update Meca to include a link between a point A and a fixed position G The force is linear: force_A = weight * ( G - A ); The point G is not associated to a Mecable, and there is no counter-force in G.

void interLink | ( | PointExact const & | pta,
|  | PointExact const & | ptb,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B The force is linear with a zero resting length:

force_A = weight * ( B - A )
force_B = weight * ( A - B )

In practice, Meca::interLink() will update the matrix mB, adding weight at the indices corresponding to A and B.

Note: with modulo, the position of the fibers may be shifted in space, and a correction is necessary to make the force calculation correct:

force_A = weight * ( B - A - offset )
force_B = weight * ( A - B + offset )

Here 'offset' is a multiple of the space periodicity, corresponding to B-A: modulo->foldOffset( A - B, offset )

In practice, Meca::interLink() will update the vector vBAS[]:

vBAS[A] += weight * offset;
vBAS[B] -= weight * offset;

In principle, what goes to vBAS[] with modulo can be derived simply by multiplying the matrix block by 'offset'.

void interLink | ( | PointInterpolated const & | pta,
|  | PointExact const & | ptb,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B The force is linear with a zero resting length: force_A = weight * ( B - A ) force_B = weight * ( A - B )

void interLink | ( | PointInterpolated const & | pta,
|  | PointInterpolated const & | ptb,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B The force is linear with a zero resting length: force_A = weight * ( B - A ) force_B = weight * ( A - B )

void interLongClamp | ( | PointExact const & | pta,
|  | Vector const & | center,
|  | real | len,
|  | real | weight
| ) |  |

Update Meca to include a non-zero resting force between A (pta) and G (center). The force is affine with non-zero resting length: force_A = (G-A) * ( len / |AG| - 1 ) There is no force on G, which is an immobile position.

void interLongLink | ( | PointExact const & | pta,
|  | PointExact const & | ptb,
|  | real | len,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B, The force is affine with non-zero resting length: force_A = weight * ( B - A ) * ( len / |AB| - 1 ) force_B = weight * ( A - B ) * ( len / |AB| - 1 )

void interLongLink | ( | PointInterpolated const & | pta,
|  | PointExact const & | pte,
|  | real | len,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B, The force is affine with non-zero resting length: force_A = weight * ( B - A ) * ( len / |AB| - 1 ) force_B = weight * ( A - B ) * ( len / |AB| - 1 )

void interLongLink | ( | PointInterpolated const & | pta,
|  | PointInterpolated const & | ptb,
|  | real | len,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B, The force is affine with non-zero resting length: force_A = weight * ( B - A ) * ( len / |AB| - 1 ) force_B = weight * ( A - B ) * ( len / |AB| - 1 )

void interPlane | ( | PointExact const & | pta,
|  | Vector const & | dir,
|  | Vector const & | g,
|  | real | weight
| ) |  |

Add interaction between pta and the plane defined by G and the normal dir The force is linear and the parallel components are removed: force(A) = weight * ( dir dir' )( G - A )

This corresponds to a frictionless plane.

The vector dir should be of norm = 1, or alternatively weight can be divided by the norm of dir.

void interPlane | ( | PointInterpolated const & | pta,
|  | Vector const & | dir,
|  | Vector const & | g,
|  | real | weight
| ) |  |

Add interaction between pta and the plane defined by G and the normal dir The force is linear and the parallel components are removed: force(A) = weight * ( dir dir' )( G - A )

This corresponds to a frictionless plane.

The vector dir should be of norm = 1, or alternatively weight can be divided by the norm of dir.

void interSideClamp | ( | PointInterpolated const & | pta,
|  | const real | g[],
|  | real | len,
|  | real | weight
| ) |  |

Update Meca to include a connection between A and a fixed position G. The force is of zero resting length, but it is taken between B and another point S which is located on the side of A: S = A + len * N, The force is linear: force_S = weight * ( G - S ) There is no counter-force in G, since G is immobile.

void interSideLink | ( | PointInterpolated const & | pta,
|  | PointExact const & | ptb,
|  | real | len,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B, but it is taken between B and a point S located on the side of A: S = A + len * N, where N is a normalized vector orthogonal to the fiber in A. S is linarly related to the two model points on the sides of an. The force is linear of zero resting length: force_S = weight * ( S - B ) force_B = weight * ( B - S )

void interSideLink | ( | PointInterpolated const & | pta,
|  | PointInterpolated const & | ptb,
|  | real | len,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B, Which is taken between B and a point S located on the side of A: S = A + len * N, where N is a normalized vector orthogonal to the fiber in an. S is linarly related to the two model points on the sides of A, P1 and P2 In 3D S is choosen in the plane of P1, P2 and B. The force is linear of zero resting length: force_S = weight * ( S - B ) force_B = weight * ( B - S )

void interSideSideLink | ( | PointInterpolated const & | pt1,
|  | PointInterpolated const & | pt2,
|  | real | len,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B, but the links are maded between SA and SB which are located on the side of A and B, respectively: SA = A + len * N_A, SB = B + len * N_B, N_X is a normalized vector orthogonal to the fiber carrying X, in X: The force is linear of zero resting length, force_SA = weight * ( SA - SB ) force_SB = weight * ( SB - SA )

void interSideSlidingLink | ( | PointInterpolated const & | pta,
|  | PointExact const & | ptb,
|  | real | len,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B, This is a combination of a SideLink with a Sliding Link: The force is linear of zero resting length, but it is taken between B, and another point S located on the side of A: S = A + len * N, where N is a normalized vector orthogonal to the fiber in A, in the direction of B. In addition, the tangential part of the force is removed.

If T is the normalized direction of the fiber in A: force_S = weight * ( 1 - T T' ) ( S - B ) force_B = weight * ( 1 - T T' ) ( B - S )

void interSideSlidingLink | ( | PointInterpolated const & | pta,
|  | PointInterpolated const & | ptb,
|  | real | len,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B, This is a combination of Side- and Sliding Links: The force is linear of zero resting length, but it is taken between B and another point S which is located on the side of A: S = A + len * N, where N is a normalized vector orthogonal to the fiber in A, in the direction of B. In addition, the part of the force tangential to A is removed.

If T is the normalized direction of the fiber in A: force_S = weight * ( 1 - T T' ) ( S - B ) force_B = weight * ( 1 - T T' ) ( B - S )

void interSlidingClamp | ( | PointInterpolated const & | pta,
|  | Vector const & | g,
|  | real | weight
| ) |  |

Update Meca to include a link between a point A and a fixed position G The force is linear and the parallel component is removed: force_A = weight * ( 1 - T T' )( G - A ) T is the vector tangent to the fiber in an. There is no counter-force in G, since G is immobile.

void interSlidingLink | ( | PointInterpolated const & | pta,
|  | PointExact const & | pte,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B, The force is linear of zero resting length, but is anisotropic: The component of the force parallel to the fiber in A is removed

If T is the normalized direction of the fiber in A: force_A = weight * ( 1 - T T' ) ( A - B ) force_B = weight * ( 1 - T T' ) ( B - A )

void interSlidingLink | ( | PointInterpolated const & | pta,
|  | PointInterpolated const & | ptb,
|  | real | weight
| ) |  |

Update Meca to include an interaction between A and B, The force is linear of zero resting length, but is anisotropic: The component of the force parallel to the fiber in A is removed

If T is the normalized direction of the fiber in A: force_A = weight * ( 1 - T T' ) ( A - B ) force_B = weight * ( 1 - T T' ) ( B - A )

void multiply | ( | const real * | X,
|  | real * | Y
| ) |  | const

calculate the matrix product needed for the conjugate gradient algorithm

Y = X - time_step * P ( mB + mC + P' ) * X;

void prepare | ( | SimulProp const * | prop | ) |

Allocate and reset matrices and vectors necessary for Meca::solve(), copy coordinates of Mecables into vPTS[]

void printMatrix | ( | std::ostream & | os | ) | const

Extract and print the full matrix associated with matVect, for debugging purposes

void solve | ( | SimulProp const * | prop,
|  | bool | precondition
| ) |  |

The equation solved is: (Xnew - Xold)/dt = P*force(X) + BrownF

Explicit integration: Xnew = Xold + dt * force(Xold) + Brown

Implicit integration using the linearized force(X) = A.X + B: ( I - dt*P*A ) ( Xnew - Xold ) = dt * P * force(Xold) + Brown

where, in both cases, Brown = sqrt(2*kT*dt*mobility) * Gaussian(0,1) Implicit integration is more robust.

Includes a constant fluid flow displacing all the objects along

## Member Data Documentation

MatrixSparseSymmetric1 mB

For interactions which have identical coefficients on the X, Y, Z subspaces

MatrixSparseSymmetric1 mC

For interactions which have different coefficients on the X, Y, Z subspaces, or which create interactions between two different subspaces.