19#ifndef CORE_OBJECT_IN_FLUID_OIF_LOCAL_FORCES_HPP
20#define CORE_OBJECT_IN_FLUID_OIF_LOCAL_FORCES_HPP
65 static constexpr int num = 3;
80 std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
85 friend boost::serialization::access;
86 template <
typename Archive>
87 void serialize(Archive &ar,
long int ) {
125 auto const dx = fp2 - fp3;
126 auto const len = dx.
norm();
127 auto const dr = len -
r0;
136 auto const lambda = len /
r0;
137 auto const spring = (std::pow(lambda, 0.5) + std::pow(lambda, -2.5)) /
138 (lambda + std::pow(lambda, -3.));
139 fac -=
ks * spring *
dr;
141 auto const f = (fac / len) * dx;
148 auto const dx = fp2 - fp3;
149 auto const len2 = dx.norm2();
150 auto const v_ij = p3.
v() - p2.
v();
176 auto const fac =
kvisc * (dx * v_ij) / len2;
177 auto const f = fac * dx;
188 auto const aa = (phi -
phi0);
190 auto const fac =
kb * aa;
199 auto const BminA = fp3 - fp2;
201 auto const factorFaNc =
202 (fp2 - fp3) * (fp1 - fp3) / BminA.norm() / Nc.norm2();
203 auto const factorFaNd =
204 (fp2 - fp3) * (fp4 - fp3) / BminA.norm() / Nd.norm2();
205 auto const factorFbNc =
206 (fp2 - fp3) * (fp2 - fp1) / BminA.norm() / Nc.norm2();
207 auto const factorFbNd =
208 (fp2 - fp3) * (fp2 - fp4) / BminA.norm() / Nd.norm2();
210 force1 -= fac * BminA.norm() / Nc.norm2() * Nc;
211 force2 += fac * (factorFaNc * Nc + factorFaNd * Nd);
212 force3 += fac * (factorFbNc * Nc + factorFbNd * Nd);
213 force4 -= fac * BminA.norm() / Nd.norm2() * Nd;
229 auto const h = (1. / 3.) * (fp1 + fp2 + fp3);
231 auto const t = sqrt(A / A0) - 1.0;
233 auto const m1 =
h - fp1;
234 auto const m2 =
h - fp2;
235 auto const m3 =
h - fp3;
237 auto const m1_length2 = m1.norm2();
238 auto const m2_length2 = m2.norm2();
239 auto const m3_length2 = m3.norm2();
242 kal * A0 * (2 * t + t * t) / (m1_length2 + m2_length2 + m3_length2);
245 force1 += (fac / 3.0) * m1;
246 force2 += (fac / 3.0) * m2;
247 force3 += (fac / 3.0) * m3;
250 handle_triangle(
kal,
A01, fp1, fp2, fp3, force1, force2, force3);
251 handle_triangle(
kal,
A02, fp2, fp3, fp4, force2, force3, force4);
253 return std::make_tuple(force2, force1, force3, force4);
Vector implementation and trait types for boost qvm interoperability.
auto unfolded_position(Utils::Vector3d const &pos, Utils::Vector3i const &image_box) const
Unfold particle coordinates to image box.
Utils::Vector< T, 3 > get_mi_vector(const Utils::Vector< T, 3 > &a, const Utils::Vector< T, 3 > &b) const
Get the minimum-image vector between two coordinates.
This file contains the defaults for ESPResSo.
#define TINY_OIF_ELASTICITY_COEFFICIENT
Tiny oif elasticity cutoff.
double angle_btw_triangles(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3, const Vector3d &P4)
Computes the angle between two triangles in 3D space.
double area_triangle(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3)
Computes the area of triangle between vectors P1,P2,P3, by computing the cross product P1P2 x P1P3 an...
Vector3d get_n_triangle(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3)
Computes the normal vector of a triangle.
Parameters for OIF local forces.
double A02
Equilibrium surface of the second triangle.
OifLocalForcesBond(double r0, double ks, double kslin, double phi0, double kb, double A01, double A02, double kal, double kvisc)
double phi0
Equilibrium angle between the two triangles.
double r0
Equilibrium bond length of triangle edges.
double kslin
Linear stretching coefficient of triangle edges.
double kb
Bending coefficient for the angle between the two triangles.
double A01
Equilibrium surface of the first triangle.
double kvisc
Viscous coefficient of the triangle vertices.
double ks
Non-linear stretching coefficient of triangle edges.
double kal
Stretching coefficient of a triangle surface.
std::tuple< Utils::Vector3d, Utils::Vector3d, Utils::Vector3d, Utils::Vector3d > calc_forces(BoxGeometry const &box_geo, Particle const &p2, Particle const &p1, Particle const &p3, Particle const &p4) const
Compute the OIF local forces.
Struct holding all information for one particle.
auto const & image_box() const