45 T G1, T G2, T v1, T v2, T scale1, T scale2, bool validationKruggelEmden) :
46 Force3D<T, PARTICLETYPE>(), _G1(G1), _G2(G2), _v1(v1), _v2(v2), _scale1(
47 scale1), _scale2(scale2), _validationKruggelEmden(validationKruggelEmden)
50 E1 = 2 * (1 + _v1) * _G1;
51 E2 = 2 * (1 + _v2) * _G2;
54 eE = (1 - util::pow(_v1, 2)) / E1 + (1 - util::pow(_v2, 2)) / E2;
58 eG = (2.0 - _v1) / _G1 + (2 - _v2) / _G2;
76 typename std::deque<PARTICLETYPE<T> >::iterator p, int pInt,
77 ParticleSystem3D<T, PARTICLETYPE>& pSys, T force[3])
80 std::vector<std::pair<size_t, T>> ret_matches;
82 pSys.getContactDetection()->getMatches(pInt, ret_matches);
84 PARTICLETYPE<T>* p2 = NULL;
87 for (
const auto& it : ret_matches) {
89 if (!util::nearZero(it.second)) {
94 T delta = (p2->getRad() + p->getRad()) - util::sqrt(it.second);
97 T M = p->getMass() * p2->getMass() / (p->getMass() + p2->getMass());
99 T R = p->getRad() * p2->getRad() / (p->getRad() + p2->getRad());
101 std::vector < T > _velR(3, T());
102 _velR[0] = -(p2->getVel()[0] - p->getVel()[0]);
103 _velR[1] = -(p2->getVel()[1] - p->getVel()[1]);
104 _velR[2] = -(p2->getVel()[2] - p->getVel()[2]);
106 std::vector < T > _d(3, T());
107 std::vector < T > _normal(3, T());
110 _d[0] = p2->getPos()[0] - p->getPos()[0];
111 _d[1] = p2->getPos()[1] - p->getPos()[1];
112 _d[2] = p2->getPos()[2] - p->getPos()[2];
114 if ( !util::nearZero(util::norm(_d)) ) {
115 _normal = util::normalize(_d);
124 T dot = velR_[0] * _normal[0] + velR_[1] * _normal[1] + velR_[2] * _normal[2];
128 std::vector < T > _velN(3, T());
129 _velN[0] = dot * _normal[0];
130 _velN[1] = dot * _normal[1];
131 _velN[2] = dot * _normal[2];
135 std::vector < T > _velT(3, T());
136 _velT[0] = _velR[0] - _velN[0];
137 _velT[1] = _velR[1] - _velN[1];
138 _velT[2] = _velR[2] - _velN[2];
145 T A =
M_PI * util::pow(util::min(p->getRad(),p2->getRad()), 2.);
146 T kn = util::pow(p->getRad()/(A*E1) + p2->getRad()/(A*E2), -1.);
150 std::vector < T > Fs_n(3, T());
151 Fs_n[0] = -kn * delta * _normal[0];
152 Fs_n[1] = -kn * delta * _normal[1];
153 Fs_n[2] = -kn * delta * _normal[2];
160 T eta_n = 0.3 * util::sqrt(4.5 * M * util::sqrt(delta) * kn);
161 if (_validationKruggelEmden) {
165 std::vector < T > Fd_n(3, T());
166 Fd_n[0] = -eta_n * _velN[0] * util::sqrt(delta);
167 Fd_n[1] = -eta_n * _velN[1] * util::sqrt(delta);
168 Fd_n[2] = -eta_n * _velN[2] * util::sqrt(delta);
170 std::vector < T > F_n(3, T());
171 F_n[0] = Fs_n[0] + Fd_n[0];
172 F_n[1] = Fs_n[1] + Fd_n[1];
173 F_n[2] = Fs_n[2] + Fd_n[2];
178 T kt = 2 * util::sqrt(2 * R) * _G1 / (2 - _v1) * util::pow(delta, 0.5);
181 T eta_t = 2 * util::sqrt(2. / 7. * M * kt);
184 std::vector < T > F_t(3, T());
185 F_t[0] = -eta_t * _velT[0];
186 F_t[1] = -eta_t * _velT[1];
187 F_t[2] = -eta_t * _velT[2];
191 force[0] = _scale1 * F_n[0] + _scale2 * F_t[0];
192 force[1] = _scale1 * F_n[1] + _scale2 * F_t[1];
193 force[2] = _scale1 * F_n[2] + _scale2 * F_t[2];
195 p->getForce()[0] += force[0] * 0.5 ;
196 p->getForce()[1] += force[1] * 0.5 ;
197 p->getForce()[2] += force[2] * 0.5 ;
198 p2->getForce()[0] -= force[0] * 0.5 ;
199 p2->getForce()[1] -= force[1] * 0.5 ;
200 p2->getForce()[2] -= force[2] * 0.5 ;