test cpu was edded...
found a bug in the oif code. It concerned the bending forces implementation. Namely, we compute the angle between two triangles in a function angle_btw_triangles (located in utils.h). The procedure was, that we found out the normals n1, n2 to the triangles, then we computed the angle between the normals and then we obtained the angle between the triangles. The step of obtaining the angle btw two normals was done from the definition of the cosine of an angle btw two vectors:

cos(a) = <n1,n2>/(|n1|.|n2|)

where <.,.> is scalar product of two vectors, and |.| is the norm of the vector. Simple expresion, however, when I computed in an usual way, that is

1. compute <n1,n2>
2. compute |n1|
3. compute |n2|
4. divide <n1,n2> by |n1| and then by |n2|

I obtained large rounding errors, around 10^-5, which, for angles close to 0, resulted in obtaining something like cos(a) = 1.000005, which is nonsense, for obvious reasons. My solution from the old days was to say simple cos(a) = 1.0.

This however led to tremendous rounding errors. After 100 time iterations, the positions of blood cell with 400 nodes differ by 0.6/400, which is unacceptable.

The remedy for this?

Strangly, I just computed it other way:
1. compute (<n1,n2>)^2
2. compute |n1|^2 - notice, here, you don't do sqrt!
3. compute |n2^2 - detto|
4. divide (<n1,n2>)^2 by |n1|^2 and then by |n2|^2
5. do sqrt of the result.

Suddenly, all erorrs went out.