192 {
193
194 auto get_tensor = [&t_w_vee](
auto a,
auto diff_a,
auto diff_diff_a,
auto b,
195 auto diff_b, auto diff_diff_b) {
201
203
206
208 auto t_hat =
getHat(t_w_vee);
209 t_diff_diff_exp(
i,
j,
k,
m) =
210
211 diff_a * FTensor::levi_civita<int>(
i,
j,
k) * t_w_vee(
m)
212
213 +
214
215 diff_a * (
216
218 FTensor::levi_civita<int>(
i,
j,
m) * t_w_vee(
k)
219
220 )
221
222 +
223
224 diff_diff_a * t_hat(
i,
j) * t_w_vee(
k) * t_w_vee(
m)
225
226 +
227
228 b * (FTensor::levi_civita<int>(
i,
l,
m) *
229 FTensor::levi_civita<int>(
l,
j,
k) +
230 FTensor::levi_civita<int>(
i,
l,
k) *
231 FTensor::levi_civita<int>(
l,
j,
m))
232
233 +
234
235 diff_b * ((t_hat(
i,
l) * FTensor::levi_civita<int>(
l,
j,
k) +
236 FTensor::levi_civita<int>(
i,
l,
k) * t_hat(
l,
j)) *
238
239 +
240
241 diff_b *
242 (
243
245
246 +
247
248 FTensor::levi_civita<int>(
i,
l,
m) * t_hat(
l,
j) * t_w_vee(
k)
249
250 +
251
252 t_hat(
i,
l) * FTensor::levi_civita<int>(
l,
j,
m) * t_w_vee(
k)
253
254 )
255
256 +
257
258 diff_diff_b * t_hat(
i,
l) * t_hat(
l,
j) * t_w_vee(
k) * t_w_vee(
m);
259
260 return t_diff_diff_exp;
261 };
262
263 if (std::fabs(theta) < std::numeric_limits<T2>::epsilon()) {
264 return get_tensor(1., -1. / 3., 1. / 15., 0., 1. / 1.2, 1. / 90.);
265 }
266
267 const auto ss = sin(theta);
268 const auto a = ss / theta;
269
270 const auto theta2 = theta * theta;
271 const auto cc = cos(theta);
272 const auto diff_a = (theta * cc - ss) / (theta2 * theta);
273 const auto diff_diff_a =
274 -(3 * theta * cc + (-3 * theta2) * ss) / pow(theta, 5);
275
276 const auto ss_2 = sin(theta / 2.);
277 const auto cc_2 = cos(theta / 2.);
278 const auto b = 2. * ss_2 * ss_2 / theta2;
279 const auto diff_b =
280 (2. * theta * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (theta2 * theta2);
281 const auto diff_diff_b =
282 (8. + (-8. + theta2) * cc - 5. * theta * ss) / (std::pow(theta, 6));
283
284 return get_tensor(
a, diff_a, diff_diff_a, b, diff_b, diff_diff_b);
285 }
#define FTENSOR_INDEX(DIM, I)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'm', 3 > m
static auto getHat(T &&w1, T &&w2, T &&w3)