81 auto const l = vec2.
norm();
84 auto const lp = vec1.
norm();
92 auto const cosPhi = (vec1 * vec2) / (lp * l);
94 auto const sinPhi = vecpro.norm() / (l * lp);
97 const double Dxx = lp /
lp0;
99 const double Dyx = 0.0;
100 const double Dyy = l /
l0 * sinPhi /
sinPhi0;
104 const double Gxy = Dxx * Dxy + Dyx * Dyy;
105 const double Gyx = Dxx * Dxy + Dyy * Dyx;
109 const double i1 = (Gxx + Gyy) - 2;
110 const double i2 = ((Gxx * Gyy) - (Gxy * Gyx)) - 1;
118 dEdI2 = -
k1 / (6.0 * (i2 + 1.0) * (i2 + 1.0));
121 dEdI1 =
k1 * (i1 + 1) / 6.0;
122 dEdI2 = -
k1 / 6.0 +
k2 * i2 / 6.0;
126 const double dI1dGxx = 1;
127 const double dI1dGxy = 0;
128 const double dI1dGyx = 0;
129 const double dI1dGyy = 1;
131 const double dI2dGxx = Gyy;
132 const double dI2dGxy = -Gyx;
134 const double dI2dGyx = -Gxy;
136 const double dI2dGyy = Gxx;
139 const double dGxxdV1x = 2 *
a1 * Dxx;
140 const double dGxxdV1y = 0;
141 const double dGxxdV2x = 2 *
a2 * Dxx;
142 const double dGxxdV2y = 0;
144 const double dGxydV1x =
a1 * Dxy +
b1 * Dxx;
145 const double dGxydV1y =
a1 * Dyy;
146 const double dGxydV2x =
a2 * Dxy +
b2 * Dxx;
147 const double dGxydV2y =
a2 * Dyy;
149 const double dGyxdV1x =
a1 * Dxy +
b1 * Dxx;
150 const double dGyxdV1y =
a1 * Dyy;
151 const double dGyxdV2x =
a2 * Dxy +
b2 * Dxx;
152 const double dGyxdV2y =
a2 * Dyy;
154 const double dGyydV1x = 2 *
b1 * Dxy;
155 const double dGyydV1y = 2 *
b1 * Dyy;
156 const double dGyydV2x = 2 *
b2 * Dxy;
157 const double dGyydV2y = 2 *
b2 * Dyy;
168 f1_rot[0] = -(dEdI1 * dI1dGxx * dGxxdV1x) - (dEdI1 * dI1dGxy * dGxydV1x) -
169 (dEdI1 * dI1dGyx * dGyxdV1x) - (dEdI1 * dI1dGyy * dGyydV1x) -
170 (dEdI2 * dI2dGxx * dGxxdV1x) - (dEdI2 * dI2dGxy * dGxydV1x) -
171 (dEdI2 * dI2dGyx * dGyxdV1x) - (dEdI2 * dI2dGyy * dGyydV1x);
172 f1_rot[1] = -(dEdI1 * dI1dGxx * dGxxdV1y) - (dEdI1 * dI1dGxy * dGxydV1y) -
173 (dEdI1 * dI1dGyx * dGyxdV1y) - (dEdI1 * dI1dGyy * dGyydV1y) -
174 (dEdI2 * dI2dGxx * dGxxdV1y) - (dEdI2 * dI2dGxy * dGxydV1y) -
175 (dEdI2 * dI2dGyx * dGyxdV1y) - (dEdI2 * dI2dGyy * dGyydV1y);
176 f2_rot[0] = -(dEdI1 * dI1dGxx * dGxxdV2x) - (dEdI1 * dI1dGxy * dGxydV2x) -
177 (dEdI1 * dI1dGyx * dGyxdV2x) - (dEdI1 * dI1dGyy * dGyydV2x) -
178 (dEdI2 * dI2dGxx * dGxxdV2x) - (dEdI2 * dI2dGxy * dGxydV2x) -
179 (dEdI2 * dI2dGyx * dGyxdV2x) - (dEdI2 * dI2dGyy * dGyydV2x);
180 f2_rot[1] = -(dEdI1 * dI1dGxx * dGxxdV2y) - (dEdI1 * dI1dGxy * dGxydV2y) -
181 (dEdI1 * dI1dGyx * dGyxdV2y) - (dEdI1 * dI1dGyy * dGyydV2y) -
182 (dEdI2 * dI2dGxx * dGxxdV2y) - (dEdI2 * dI2dGxy * dGxydV2y) -
183 (dEdI2 * dI2dGyx * dGyxdV2y) - (dEdI2 * dI2dGyy * dGyydV2y);
190 auto forces = RotateForces(f1_rot, f2_rot, vec1, vec2);
std::tuple< Utils::Vector3d, Utils::Vector3d, Utils::Vector3d > RotateForces(Utils::Vector2d const &f1_rot, Utils::Vector2d const &f2_rot, Utils::Vector3d const &v12, Utils::Vector3d const &v13)
Rotate calculated trielastic forces in the 2d plane back to the 3d plane.
IBMTriel(int ind1, int ind2, int ind3, double maxDist, tElasticLaw elasticLaw, double k1, double k2)
Set the IBM Triel parameters.