p2p_cylinder_collisions.hpp 7.76 KB
Newer Older
1
2
3
4
/******************************************************************************
*                                                                             *
*  Copyright 2019 Max Planck Institute for Dynamics and Self-Organization     *
*                                                                             *
Cristian Lalescu's avatar
Cristian Lalescu committed
5
*  This file is part of TurTLE.                                               *
6
*                                                                             *
Cristian Lalescu's avatar
Cristian Lalescu committed
7
*  TurTLE is free software: you can redistribute it and/or modify             *
8
9
10
11
*  it under the terms of the GNU General Public License as published          *
*  by the Free Software Foundation, either version 3 of the License,          *
*  or (at your option) any later version.                                     *
*                                                                             *
Cristian Lalescu's avatar
Cristian Lalescu committed
12
*  TurTLE is distributed in the hope that it will be useful,                  *
13
14
15
16
17
*  but WITHOUT ANY WARRANTY; without even the implied warranty of             *
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
*  GNU General Public License for more details.                               *
*                                                                             *
*  You should have received a copy of the GNU General Public License          *
Cristian Lalescu's avatar
Cristian Lalescu committed
18
*  along with TurTLE.  If not, see <http://www.gnu.org/licenses/>             *
19
20
21
22
*                                                                             *
* Contact: Cristian.Lalescu@ds.mpg.de                                         *
*                                                                             *
******************************************************************************/
23
24
#ifndef P2P_CYLINDER_GHOST_COLLISIONS_HPP
#define P2P_CYLINDER_GHOST_COLLISIONS_HPP
25
26

#include <cstring>
27
#include <math.h>       /* for sqrt, abs */
28

29
30
31
#include <set>
#include <utility>
#include <vector>
32
#include "particles/p2p_ghost_collision_base.hpp"
33
34


35
template <class real_number, class partsize_t>
36
37
38
39
40
41
class p2p_cylinder_ghost_collisions: public p2p_ghost_collision_base<real_number, partsize_t>
{
    private:
        // Following doubles are needed for the collision computation
        double width;
        double length;
42
43

public:
44
    p2p_cylinder_ghost_collisions() : width(1.),length(1.){}
45
46

    // Copy constructor use a counter set to zero
47
    p2p_cylinder_ghost_collisions(const p2p_cylinder_ghost_collisions&){}
48
49

    template <int size_particle_positions, int size_particle_rhs>
50
    void compute_interaction(const partsize_t idx_part1,
51
                             const real_number pos_part1[], real_number /*rhs_part1*/[],
52
                             const partsize_t idx_part2,
53
                             const real_number pos_part2[], real_number /*rhs_part2*/[],
54
                             const real_number /*dist_pow2*/,  const real_number /*cutoff*/,
55
                             const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/){
56
        double x, y, z, pq, xp, xq, t, s, x_dist, y_dist, z_dist ,min_distance, pi, pi2;
57

58
59
        pi2 = atan(1.0)*8.0;
        pi = atan(1.0)*4.0;
60
61

        /* Relative position vector of the two particles (x,y,z)^T */
Cristian Lalescu's avatar
Cristian Lalescu committed
62
        /* complicated usage of fmod, fabs and sgn because C++ fmod does not do what it should. */
63
64
65
        x = std::fmod(pos_part2[0],pi2)-std::fmod(pos_part1[0],pi2);
        y = std::fmod(pos_part2[1],pi2)-std::fmod(pos_part1[1],pi2);
        z = std::fmod(pos_part2[2],pi2)-std::fmod(pos_part1[2],pi2);
66

67
68
69
        x = ( std::fmod( std::fabs(x) +pi ,pi2) - pi ) * sgn(x) ;
        y = ( std::fmod( std::fabs(y) +pi ,pi2) - pi ) * sgn(y) ;
        z = ( std::fmod( std::fabs(z) +pi ,pi2) - pi ) * sgn(z) ;
Cristian Lalescu's avatar
Cristian Lalescu committed
70

71
72
        if( sqrt(x*x+y*y+z*z) <= length )
        {
Cristian Lalescu's avatar
Cristian Lalescu committed
73

74
75
76
77
78
79
        /* p and q are the orientation vectors of the first and second particles. */
        /* pq, xp, xq are scalar products of these vectors with x, relative position */
        pq = pos_part1[3] * pos_part2[3] + pos_part1[4] * pos_part2[4] + pos_part1[5] * pos_part2[5];
        xp = x * pos_part1[3] + y * pos_part1[4] + z * pos_part1[5];
        xq = x * pos_part2[3] + y * pos_part2[4] + z * pos_part2[5];
        /* t and s parametrize the two rods. Find min distance: */
80
81
82
        assert(this->length > 0);
        t = 2.0/(this->length*(pq*pq-1.0))*(-xp+pq*xq);
        s = 2.0/(this->length*(pq*pq-1.0))*(-pq*xp+xq);
83
84
85
        /* Test if -1<s<1 and -1<t<1 */
        if( abs(t)<=1.0 and abs(s)<=1.0 ) {
            /* Get minimal distance in case of both t and s in {-1,1}. Else: check edges */
86
87
88
            x_dist = this->length*0.5*t*pos_part1[3]-this->length*0.5*s*pos_part2[3]-x;
            y_dist = this->length*0.5*t*pos_part1[4]-this->length*0.5*s*pos_part2[4]-y;
            z_dist = this->length*0.5*t*pos_part1[5]-this->length*0.5*s*pos_part2[5]-z;
89
            min_distance = sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist);
90
91
92
        }
        else
        {
93
            min_distance = this->length;
94
95
            /* t fixed at 1, find min along s */
            t = 1.0;
96
            s = t*pq-2.0/this->length*xq;
97
            if( abs(s)>1.0 ) { s = s / abs(s) ;}
98
99
100
            x_dist = this->length*0.5*t*pos_part1[3]-this->length*0.5*s*pos_part2[3]-x;
            y_dist = this->length*0.5*t*pos_part1[4]-this->length*0.5*s*pos_part2[4]-y;
            z_dist = this->length*0.5*t*pos_part1[5]-this->length*0.5*s*pos_part2[5]-z;
101
102
103
            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
            /* t fixed at -1, find min along s */
            t = -1.0;
104
            s = t*pq-2.0/this->length*xq;
105
            if( abs(s)>1.0 ) { s = s / abs(s) ;}
106
107
108
            x_dist = this->length*0.5*t*pos_part1[3]-this->length*0.5*s*pos_part2[3]-x;
            y_dist = this->length*0.5*t*pos_part1[4]-this->length*0.5*s*pos_part2[4]-y;
            z_dist = this->length*0.5*t*pos_part1[5]-this->length*0.5*s*pos_part2[5]-z;
109
110
111
            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
            /* s fixed at 1, find min along t */
            s = 1.0;
112
            t = s*pq+2.0/this->length*xp;
113
            if( abs(t)>1.0 ) { t = t / abs(t) ;}
114
115
116
            x_dist = this->length*0.5*t*pos_part1[3]-this->length*0.5*s*pos_part2[3]-x;
            y_dist = this->length*0.5*t*pos_part1[4]-this->length*0.5*s*pos_part2[4]-y;
            z_dist = this->length*0.5*t*pos_part1[5]-this->length*0.5*s*pos_part2[5]-z;
117
118
119
            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
            /* s fixed at -1, find min along t */
            s = -1.0;
120
            t = s*pq+2.0/this->length*xp;
121
            if( abs(t)>1.0 ) { t = t / abs(t) ;}
122
123
124
            x_dist = this->length*0.5*t*pos_part1[3]-this->length*0.5*s*pos_part2[3]-x;
            y_dist = this->length*0.5*t*pos_part1[4]-this->length*0.5*s*pos_part2[4]-y;
            z_dist = this->length*0.5*t*pos_part1[5]-this->length*0.5*s*pos_part2[5]-z;
125
126
127
            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
        }
        /* If cylinders overlap count it as a collision */
Tobias Baetge's avatar
Tobias Baetge committed
128
129
130
131
132
133
        if( min_distance<=width ){
            std::pair <partsize_t, partsize_t> single_collision_pair(idx_part1, idx_part2);
            this->collision_pairs_local.insert(single_collision_pair);
            //DEBUG_MSG("inside compute interaction idx_part1 = %ld and idx_part2 = %ld\n", idx_part1, idx_part2);
            assert(idx_part1!=idx_part2);
        }
134

135
		}
136
137
    }

138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
    void set_width(const double WIDTH)
    {
        this->width = WIDTH;
    }

    double get_width()
    {
        return this->width;
    }

    void set_length(const double LENGTH)
    {
        this->length = LENGTH;
    }

    double get_length()
    {
        return this->length;
    }
157
158
159
};


160
#endif // P2P_CYLINDER_GHOST_COLLISIONS_HPP
161