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
212
213 +
214
215 diff_a * (
216
219
220 )
221
222 +
223
224 diff_diff_a * t_hat(
i,
j) * t_w_vee(
k) * t_w_vee(
m)
225
226 +
227
232
233 +
234
238
239 +
240
241 diff_b *
242 (
243
245
246 +
247
249
250 +
251
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
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
FTensor::Index< 'm', 3 > m
static auto getHat(T &&w1, T &&w2, T &&w3)