Oif corrections

[RudolfWeeber]

This ports back fixes from the object in fluid development branch. In particular, the bending force between two triangeles becomes torque-free with this.
The test in this PR fails without the correctiosn in oif_local_force.hpp

[RudolfWeeber]

@icimrak could you please re-check this. There were some merge conflicts I had to resolve.

[codecov]

Codecov Report

Merging #3385 into python will decrease coverage by <1%.
The diff coverage is 84%.

@@          Coverage Diff           @@
##           python   #3385   +/-   ##
======================================
- Coverage      87%     87%   -1%     
======================================
  Files         536     534    -2     
  Lines       24403   24356   -47     
======================================
- Hits        21252   21202   -50     
- Misses       3151    3154    +3

Impacted Files	Coverage Δ
src/core/io/writer/h5md_core.cpp	89% <ø> (ø)	⬆️
src/core/nsquare.cpp	100% <ø> (ø)	⬆️
src/core/grid_based_algorithms/lbgpu.cpp	91% <ø> (ø)	⬆️
src/core/layered.cpp	79% <ø> (ø)	⬆️
src/core/communication.cpp	94% <ø> (ø)	⬆️
src/core/event.cpp	94% <ø> (ø)	⬆️
src/core/virtual_sites.cpp	85% <ø> (ø)	⬆️
src/core/reaction_ensemble.hpp	78% <ø> (ø)	⬆️
src/core/immersed_boundary/ImmersedBoundaries.cpp	20% <0%> (ø)	⬆️
src/core/io/mpiio/mpiio.cpp	84% <100%> (-1%)	⬇️
... and 23 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update f7726af...cfbb876. Read the comment docs.

[icimrak]

These changes are ok....

[RudolfWeeber]

@jngrad, I chrrey-picked your suggestion. I also reverted the randomization of the test, since it is not easy to guarantee that this doesn't flip some of the normal vectors of the triangles making up the surface. there is a deterministic anisotropic deformation now, instead.

As for the corner cases: It is not clear to me how to construct triangles such that the dot product of the normal vectors is >1 or <-1. Formally, that never happens. One would have to construct triangles with normal vectors triggering a rounding error after normalization.

@icimrak, could you suggest a docstring explaining the convention for the resulting angles, and supply two sets of two triangles with angles in the two half-spaces? Please avoid zeros and ones in the coordinates.
I'll then write down the unit test.

[fweik]

@RudolfWeeber the dot product is (exactly) +/- 1 if both triangles have the same points, in which case the normals are either parallel or anti-parallel and hence their inner product is +/- 1.

[RudolfWeeber]

@fweik, my understanding of https://en.cppreference.com/w/cpp/numeric/math/acos is that +,- 1.0 are permissible. Due to floating point rounding, ocassionally 1 + epsilon occurs, which is not permissible and has to be capped. So, if we want to test the corner cases, that would be them.

[fweik]

You were saying that +/-1 are formally impossible, which is not correct. Rather, this corresponds to the very common case where both triangles are in the same plane.

[fweik]

LGTM, would be perfect if there was a test for the revised angle_btw_triangles function.

[icimrak]

Hi Rudolf,

Do we get in touch via skype? We can then discuss issues listed on github, as well as the article….

Here are the times we are available:

9.1. only between 10am till 12pm
10.1. between 1pm and 3 30pm

Please let us know when are you available,

Thanks,
Ivan

[icimrak]

Here is the suggestion for a docstring for function call:
angle_btw_triangles(P1,P2,P3,P4)

Returns the angle between two triangles in 3D space given by points P1P2P3 and P2P3P4. Note, that the common edge is given as the second and the third argument. Here, the angle can have values from 0 to 2 * PI, depending on the orientation of the two triangles. So the angle can be convex or concave. For each triangle, an inward direction has been defined as the direction of one of the two normal vectors. Particularly, for triangle P1P2P3 it is the vector N1 = P2P1 x P2P3 and for triangle P2P3P4 it is N2 = P2P3 x P2P4. The method first computes the angle between N1 and N2, which gives always value between 0 and PI and then it checks whether this value must be corrected to a value between PI and 2 * PI.

As an example, consider 4 points A,B,C,D in space given by coordinates A = [1,1,1], B = [2,1,1], C = [1,2,1], D = [1,1,2]. We want to determine the angle between triangles ABC and ACD. In case that the orientations of the triangle ABC is [0,0,1] and orientation of ACD is [1,0,0], then the resulting angle must be PI/2.0. To get correct result, note that the common edge is AC, and one must call the method as
angle_btw_triangles(B,A,C,D).
With this call we have ensured that N1 = AB x AC (which coincides with [0,0,1]) and N2 = AC x AD (which coincides with [1,0,0]).

Alternatively, if the orientations of the two triangles were the oppisite, the correct call would be angle_btw_triangles(B,C,A,D)
so that N1 = CB x CA and N2 = CA x CD.

[fweik]

Can be merged if

 /builds/espressomd/espresso/src/utils/tests/triangle_functions_test.cpp:90:2: error: extra ';' [-Werror=pedantic]

is fixed.