16  typedef typename std::remove_pointer<T>::type 
Type;
 
 
   19  typedef typename std::remove_pointer<T>::type 
Type;
 
 
   27  template <
typename T> 
inline static auto getVee(T &&w1, T &&w2, T &&w3) {
 
   30        std::forward<T>(w1), std::forward<T>(w2), std::forward<T>(w3)
 
 
   35  template <
typename T> 
inline static auto getHat(T &&w1, T &&w2, T &&w3) {
 
 
   39  template <
typename A> 
inline static auto getVee(
A &&t_w_hat) {
 
 
   43  template <
typename A> 
inline static auto getHat(
A &&t_w_vee) {
 
 
   47  template <
typename A, 
typename B>
 
   48  inline static auto exp(
A &&t_w_vee, 
B &&theta) {
 
   49    return expImpl(std::forward<A>(t_w_vee), std::forward<B>(theta));
 
 
   52  template <
typename A, 
typename B>
 
   53  inline static auto Jl(
A &&t_w_vee, 
B &&theta) {
 
   54    return JlImpl(std::forward<A>(t_w_vee), std::forward<B>(theta));
 
 
   57  template <
typename A, 
typename B>
 
   58  inline static auto Jr(
A &&t_w_vee, 
B &&theta) {
 
   59    return JrImpl(std::forward<A>(t_w_vee), std::forward<B>(theta));
 
 
   62  template <
typename A, 
typename B, 
typename C>
 
   63  inline static auto action(
A &&t_w_vee, 
B &&theta, C &&t_A) {
 
   64    return actionImpl(std::forward<A>(t_w_vee), std::forward<B>(theta),
 
   65                      std::forward<C>(t_A));
 
 
   68  template <
typename A, 
typename B>
 
   70    return dActionJlImpl(std::forward<A>(t_w_vee), std::forward<B>(t_A));
 
 
   73  template <
typename A, 
typename B>
 
   75    return dActionJrImpl(std::forward<A>(t_w_vee), std::forward<B>(t_A));
 
 
   78  template <
typename A, 
typename B>
 
   79  inline static auto diffExp(
A &&t_w_vee, 
B &&theta) {
 
   80    return diffExpImpl(std::forward<A>(t_w_vee), std::forward<B>(theta));
 
 
   83  template <
typename A, 
typename B>
 
   85    return diffDiffExpImpl(std::forward<A>(t_w_vee), std::forward<B>(theta));
 
 
   89  inline static constexpr int dim = 3;
 
  101    t_w_vee(
k) = (levi_civita(
i, 
j, 
k) * t_w_hat(
i, 
j)) / 2;
 
 
  105  template <
typename T>
 
  109    t_w_hat(
i, 
j) = levi_civita(
i, 
j, 
k) * t_w_vee(
k);
 
 
  113  template <
typename T1, 
typename T2, 
typename T3>
 
  115                                     const T2 alpha, 
const T3 beta) {
 
  118    auto t_hat = 
getHat(t_w_vee);
 
  120                beta * (t_hat(
i, 
k) * t_hat(
k, 
j));
 
 
  124  template <
typename T1, 
typename T2>
 
  127    if (std::fabs(theta) < std::numeric_limits<T2>::epsilon()) {
 
  130    const auto s = sin(theta);
 
  131    const auto s_half = sin(theta / 2);
 
  132    const auto a = s / theta;
 
  133    const auto b = 2 * (s_half / theta) * (s_half / theta);
 
 
  137  template <
typename T1, 
typename T2>
 
  141    auto get_tensor = [&t_w_vee](
auto a, 
auto diff_a, 
auto b, 
auto diff_b) {
 
  149      auto t_hat = 
getHat(t_w_vee);
 
  150      t_diff_exp(
i, 
j, 
k) =
 
  152          a * FTensor::levi_civita<int>(
i, 
j, 
k)
 
  156          diff_a * t_hat(
i, 
j) * t_w_vee(
k)
 
  160          b * (t_hat(
i, 
l) * FTensor::levi_civita<int>(
l, 
j, 
k) +
 
  161               FTensor::levi_civita<int>(
i, 
l, 
k) * t_hat(
l, 
j))
 
  165          diff_b * t_hat(
i, 
l) * t_hat(
l, 
j) * t_w_vee(
k);
 
  170    if(std::fabs(theta) < std::numeric_limits<T2>::epsilon()){
 
  171      return get_tensor(1., -1. / 3., 0., 1. / 1.2);
 
  174    const auto ss = sin(theta);
 
  175    const auto a = ss / theta;
 
  177    const auto theta2 = theta * theta;
 
  178    const auto cc = cos(theta);
 
  179    const auto diff_a = (theta * cc - ss) / (theta2 * theta);
 
  181    const auto ss_2 = sin(theta / 2.);
 
  182    const auto cc_2 = cos(theta / 2.);
 
  183    const auto b = 2. * ss_2 * ss_2 / theta2;
 
  185        (2. * theta * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (theta2 * theta2);
 
  187    return get_tensor(
a, diff_a, b, diff_b);
 
 
  190  template <
typename T1, 
typename T2>
 
  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) {
 
  208      auto t_hat = 
getHat(t_w_vee);
 
  209      t_diff_diff_exp(
i, 
j, 
k, 
m) =
 
  211          diff_a * FTensor::levi_civita<int>(
i, 
j, 
k) * t_w_vee(
m)
 
  218                       FTensor::levi_civita<int>(
i, 
j, 
m) * t_w_vee(
k)
 
  224          diff_diff_a * t_hat(
i, 
j) * t_w_vee(
k) * t_w_vee(
m)
 
  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))
 
  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)) *
 
  248                  FTensor::levi_civita<int>(
i, 
l, 
m) * t_hat(
l, 
j) * t_w_vee(
k)
 
  252                  t_hat(
i, 
l) * FTensor::levi_civita<int>(
l, 
j, 
m) * t_w_vee(
k)
 
  258          diff_diff_b * t_hat(
i, 
l) * t_hat(
l, 
j) * t_w_vee(
k) * t_w_vee(
m);
 
  260      return t_diff_diff_exp;
 
  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.);
 
  267    const auto ss = sin(theta);
 
  268    const auto a = ss / theta;
 
  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);
 
  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;
 
  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));
 
  284    return get_tensor(
a, diff_a, diff_diff_a, b, diff_b, diff_diff_b);
 
 
  287  template <
typename T1, 
typename T2>
 
  290    if (std::fabs(theta) < std::numeric_limits<T2>::epsilon()) {
 
  293    const auto s = sin(theta);
 
  294    const auto s_half = sin(theta / 2);
 
  295    const auto a = 2 * (s_half / theta) * (s_half / theta);
 
  296    const auto b = ((theta - s) / theta) / theta / theta;
 
 
  300  template <
typename T1, 
typename T2>
 
  303    if (std::fabs(theta) < std::numeric_limits<T2>::epsilon()) {
 
  306    const auto s = sin(theta);
 
  307    const auto s_half = sin(theta / 2);
 
  308    const auto a = 2 * (s_half / theta) * (s_half / theta);
 
  309    const auto b = ((theta - s) / theta) / theta / theta;
 
 
  313  template <
typename T1, 
typename T2, 
typename T3>
 
  319    t_B(
i, 
j) = 
exp(t_w_vee, theta)(
i,
k) * t_A(
k, 
j);
 
 
 
 
#define FTENSOR_INDEX(DIM, I)
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
Tensors class implemented by Walter Landry.
FTensor::Index< 'm', 3 > m
static auto diffDiffExp(A &&t_w_vee, B &&theta)
static auto diffExpImpl(FTensor::Tensor1< T1, dim > &t_w_vee, const T2 theta)
static auto dActionJr(A &&t_w_vee, B &&t_A)
static auto Jr(A &&t_w_vee, B &&theta)
static auto Jl(A &&t_w_vee, B &&theta)
static auto action(A &&t_w_vee, B &&theta, C &&t_A)
static auto JrImpl(FTensor::Tensor1< T1, dim > &t_w_vee, const T2 theta)
static FTENSOR_INDEX(dim, l)
static FTENSOR_INDEX(dim, k)
static auto JlImpl(FTensor::Tensor1< T1, dim > &t_w_vee, const T2 &theta)
static auto getVee(A &&t_w_hat)
static auto getHatImpl(FTensor::Tensor1< T, dim > &t_w_vee)
static auto getHat(T &&w1, T &&w2, T &&w3)
static FTENSOR_INDEX(dim, j)
static auto genericFormImpl(FTensor::Tensor1< T1, dim > &t_w_vee, const T2 alpha, const T3 beta)
static auto diffExp(A &&t_w_vee, B &&theta)
static auto diffDiffExpImpl(FTensor::Tensor1< T1, dim > &t_w_vee, const T2 theta)
static auto dActionJl(A &&t_w_vee, B &&t_A)
static auto actionImpl(FTensor::Tensor1< T1, dim > &t_w_vee, const T2 theta, FTensor::Tensor2_symmetric< T3, dim > &t_A)
static auto expImpl(FTensor::Tensor1< T1, dim > &t_w_vee, const T2 theta)
static auto getHat(A &&t_w_vee)
static FTENSOR_INDEX(dim, n)
static FTENSOR_INDEX(dim, m)
static auto getVeeImpl(FTensor::Tensor2< T, dim, dim > &t_w_hat)
static FTENSOR_INDEX(dim, i)
static auto getVee(T &&w1, T &&w2, T &&w3)
static auto exp(A &&t_w_vee, B &&theta)