#include "src/ftensor/src/MatrixFunctionTemplate.hpp"

Collaboration diagram for EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >:

Public Types
using	Val = const FTensor::Tensor1< T1, Dim >

using	Vec = const FTensor::Tensor2< T2, Dim, Dim >

using	Fun = boost::function< double(const double)>

using	V = double

template<int N>
using	Number = FTensor::Number< N >

template<char c>
using	I = typename FTensor::Index< c, Dim >

using	NumberNb = Number< NB >

using	NumberDim = Number< Dim >

Public Member Functions
	EigenMatrixImp (Val &t_val, Vec &t_vec)

auto	getMat (Fun f)
	Get matrix.

auto	getDiffMat (Fun f, Fun d_f)
	Get derivative of matrix.

template<typename T >
auto	getDiffDiffMat (Fun f, Fun d_f, Fun dd_f, T &t_S)
	Get second directive of matrix.

Private Attributes
Val &	tVal

Vec &	tVec

FTensor::Tensor2_symmetric< V, Dim >	aM [Dim]

FTensor::Ddg< V, Dim, Dim >	aMM [Dim][Dim]

FTensor::Ddg< V, Dim, Dim >	aG [Dim][Dim]

FTensor::Ddg< V, Dim, Dim >	aS [(Dim *(Dim+1))/2]

FTensor::Ddg< V, Dim, Dim >	aSM [(Dim - 1) Dim][(Dim (Dim+1))/2]

FTensor::Ddg< V, Dim, Dim >	d2MType0 [Dim][(Dim *(Dim+1))/2]

FTensor::Ddg< V, Dim, Dim >	d2MType1 [Dim][(Dim *(Dim+1))/2]

FTensor::Tensor2_symmetric< V, Dim >	aF2

FTensor::Tensor2< V, Dim, Dim >	aF

FTensor::Tensor1< V, Dim >	fVal

FTensor::Tensor1< V, Dim >	dfVal

FTensor::Tensor1< V, Dim >	ddfVal

Friends
template<typename E , typename C >
struct	d2MCoefficients

template<typename E , typename C , typename G >
struct	d2MImpl

template<typename E , typename C >
struct	Fdd4MImpl

template<typename E , typename C >
struct	ReconstructMatImpl

template<typename E , typename C >
struct	FirstMatrixDirectiveImpl

template<typename E , typename C >
struct	SecondMatrixDirectiveImpl

template<typename E , typename C , typename T >
struct	GetDiffMatImpl

template<typename E , typename C , typename T3 , typename T4 >
struct	GetDiffDiffMatImpl

Detailed Description

template<typename T1, typename T2, int NB, int Dim>
struct EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >

Definition at line 931 of file MatrixFunctionTemplate.hpp.

Member Typedef Documentation

◆ Fun

template<typename T1 , typename T2 , int NB, int Dim>

using EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::Fun = boost::function<double(const double)>

Definition at line 935 of file MatrixFunctionTemplate.hpp.

◆ I

template<typename T1 , typename T2 , int NB, int Dim>

template<char c>

using EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::I = typename FTensor::Index<c, Dim>

Definition at line 939 of file MatrixFunctionTemplate.hpp.

◆ Number

template<typename T1 , typename T2 , int NB, int Dim>

template<int N>

using EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::Number = FTensor::Number<N>

Definition at line 938 of file MatrixFunctionTemplate.hpp.

◆ NumberDim

template<typename T1 , typename T2 , int NB, int Dim>

using EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::NumberDim = Number<Dim>

Definition at line 942 of file MatrixFunctionTemplate.hpp.

◆ NumberNb

template<typename T1 , typename T2 , int NB, int Dim>

using EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::NumberNb = Number<NB>

Definition at line 941 of file MatrixFunctionTemplate.hpp.

◆ V

template<typename T1 , typename T2 , int NB, int Dim>

using EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::V = double

Definition at line 936 of file MatrixFunctionTemplate.hpp.

◆ Val

template<typename T1 , typename T2 , int NB, int Dim>

using EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::Val = const FTensor::Tensor1<T1, Dim>

Definition at line 933 of file MatrixFunctionTemplate.hpp.

◆ Vec

template<typename T1 , typename T2 , int NB, int Dim>

using EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::Vec = const FTensor::Tensor2<T2, Dim, Dim>

Definition at line 934 of file MatrixFunctionTemplate.hpp.

Constructor & Destructor Documentation

◆ EigenMatrixImp()

template<typename T1 , typename T2 , int NB, int Dim>

EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::EigenMatrixImp	(	Val &	t_val,
		Vec &	t_vec
	)

inline

Definition at line 944 of file MatrixFunctionTemplate.hpp.

                                         : tVal(t_val), tVec(t_vec) {
 
    FTensor::Index<'i', Dim> i;
    FTensor::Index<'j', Dim> j;
    FTensor::Index<'k', Dim> k;
    FTensor::Index<'l', Dim> l;
 
    for (auto aa = 0; aa != Dim; ++aa) {
      auto &M = aM[aa];
      for (auto ii = 0; ii != Dim; ++ii)
        for (auto jj = 0; jj <= ii; ++jj)
          M(ii, jj) = tVec(aa, ii) * tVec(aa, jj);
    }
 
    for (auto aa = 0; aa != Dim; ++aa) {
      for (auto bb = 0; bb != Dim; ++bb) {
        auto &Ma = aM[aa];
        auto &Mb = aM[bb];
        auto &MM = aMM[aa][bb];
        MM(i, j, k, l) = Ma(i, j) * Mb(k, l);
      }
    }
 
    for (auto aa = 0; aa != Dim; ++aa) {
      for (auto bb = 0; bb != Dim; ++bb) {
        if (aa != bb) {
          auto &MM = aMM[aa][bb];
          auto &G = aG[aa][bb];
          G(i, j, k, l) = MM(i, k, j, l) || MM(i, l, j, k);
        }
      }
    }
 
    for (auto aa = 0; aa != Dim; ++aa) {
      for (auto bb = 0; bb != Dim; ++bb) {
        if (aa != bb) {
          auto &Gab = aG[aa][bb];
          auto &Gba = aG[bb][aa];
          auto &S = aS[get_sym_index(aa, bb, Dim)];
          S(i, j, k, l) = Gab(i, j, k, l) + Gba(i, j, k, l);
        }
      }
    }
  }

Member Function Documentation

◆ getDiffDiffMat()

template<typename T1 , typename T2 , int NB, int Dim>

template<typename T >

auto EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::getDiffDiffMat	(	Fun	f,
		Fun	d_f,
		Fun	dd_f,
		T &	t_S
	)

inline

Get second directive of matrix.

\[ LS_{klmn} = S_{ij} \frac{\partial^2 B_{ij}}{\partial A_{kl} \partial A_{mn} } \]

Template Parameters

T

Parameters

t_val	eigen values vector
t_vec	eigen vectors matrix
f	function
d_f	derivative of function
dd_f	second derivative of function
t_S	second rank tensor S

Returns: auto second derivatives, forth order tensor with minor symmetries

Definition at line 1072 of file MatrixFunctionTemplate.hpp.

                                                               {
 
    for (auto aa = 0; aa != Dim; ++aa)
      fVal(aa) = f(tVal(aa));
 
    for (auto aa = 0; aa != Dim; ++aa)
      dfVal(aa) = d_f(tVal(aa));
 
    for (auto aa = 0; aa != Dim; ++aa)
      ddfVal(aa) = dd_f(tVal(aa));
 
    for (auto aa = 0; aa != Dim; ++aa)
      for (auto bb = 0; bb != aa; ++bb) {
        aF(aa, bb) = 1 / (tVal(aa) - tVal(bb));
        aF(bb, aa) = -aF(aa, bb);
        aF2(aa, bb) = aF(aa, bb) * aF(aa, bb);
      }
 
    FTensor::Index<'i', Dim> i;
    FTensor::Index<'j', Dim> j;
    FTensor::Index<'k', Dim> k;
    FTensor::Index<'l', Dim> l;
 
    for (auto aa = 0; aa != Dim; ++aa) {
      for (auto bb = 0; bb != Dim; ++bb) {
        if (aa != bb) {
          const auto &S = aS[get_sym_index(aa, bb, Dim)];
          const auto &M = aM[aa];
          auto &SMmn = aSM[get_nodiag_index(aa, bb, Dim)];
          for (auto mm = 0; mm != Dim; ++mm) {
            for (auto nn = mm; nn != Dim; ++nn) {
              SMmn[get_sym_index(mm, nn, Dim)](i, j, k, l) =
                  S(i, j, k, l) * M(mm, nn);
            }
          }
        }
      }
    }
 
    for (auto aa = 0; aa != Dim; ++aa) {
      for (auto mm = 0; mm != Dim; ++mm) {
        for (auto nn = mm; nn != Dim; ++nn) {
          d2MType0[aa][get_sym_index(mm, nn, Dim)](i, j, k, l) = 0;
        }
      }
    }
 
    if constexpr (NB == 3)
      for (auto aa = 0; aa != Dim; ++aa) {
        for (auto bb = 0; bb != Dim; ++bb) {
          if (aa != bb) {
            const V v = dfVal(aa) * aF(aa, bb);
            for (auto mm = 0; mm != Dim; ++mm) {
              for (auto nn = mm; nn != Dim; ++nn) {
                d2MType0[aa][get_sym_index(mm, nn, Dim)](i, j, k, l) +=
                    v * aSM[get_nodiag_index(aa, bb, Dim)]
                           [get_sym_index(mm, nn, Dim)](i, j, k, l);
              }
            }
          }
        }
      }
 
    if constexpr (NB == 2)
      for (auto aa = 0; aa != Dim; ++aa) {
        for (auto bb = 0; bb != Dim; ++bb) {
          if (aa != bb) {
            V v;
            if (aa == 1 || bb == 1)
              v = dfVal(aa) * aF(aa, bb);
            else
              v = ddfVal(aa) / 2;
            for (auto mm = 0; mm != Dim; ++mm) {
              for (auto nn = mm; nn != Dim; ++nn) {
                d2MType0[aa][get_sym_index(mm, nn, Dim)](i, j, k, l) +=
                    v * aSM[get_nodiag_index(aa, bb, Dim)]
                           [get_sym_index(mm, nn, Dim)](i, j, k, l);
              }
            }
          }
        }
      }
 
    if constexpr (NB == 1)
      for (auto aa = 0; aa != Dim; ++aa) {
        for (auto bb = 0; bb != Dim; ++bb) {
          if (aa != bb) {
            const V v = ddfVal(aa) / 2;
            for (auto mm = 0; mm != Dim; ++mm) {
              for (auto nn = mm; nn != Dim; ++nn) {
                d2MType0[aa][get_sym_index(mm, nn, Dim)](i, j, k, l) +=
                    v * aSM[get_nodiag_index(aa, bb, Dim)]
                           [get_sym_index(mm, nn, Dim)](i, j, k, l);
              }
            }
          }
        }
      }
 
    for (auto aa = 0; aa != Dim; ++aa) {
      for (auto mm = 0; mm != Dim; ++mm) {
        for (auto nn = 0; nn != Dim; ++nn) {
          if (nn != mm)
            d2MType1[mm][get_sym_index(aa, nn, Dim)](i, j, k, l) = 0;
        }
      }
    }
 
    if constexpr (NB == 3)
      for (auto aa = 0; aa != Dim; ++aa) {
        for (auto bb = 0; bb != Dim; ++bb) {
          if (aa != bb) {
            const auto &S = aS[get_sym_index(aa, bb, Dim)];
            const auto v0 = aF(aa, bb);
            for (auto cc = 0; cc != Dim; ++cc) {
              for (auto dd = 0; dd != Dim; ++dd) {
                if (cc != dd) {
                  const double v1 = fVal(cc) * aF(cc, dd);
                  d2MType1[dd][get_sym_index(aa, cc, Dim)](i, j, k, l) +=
                      (v1 * v0) * S(i, j, k, l);
                }
              }
            }
          }
        }
      }
 
    if constexpr (NB == 2)
      for (auto aa = 0; aa != Dim; ++aa) {
        for (auto bb = 0; bb != Dim; ++bb) {
          if (aa != bb)
            for (auto cc = 0; cc != Dim; ++cc) {
              for (auto dd = 0; dd != Dim; ++dd)
                if (cc != dd) {
 
                  V r;
 
                  if ((cc == 1 || dd == 1) && (aa == 1 || bb == 1))
                    r = fVal(cc) * aF(cc, dd) * aF(aa, bb);
                  else if (cc != 1 && dd != 1 && aa != 1 && bb != 1) {
 
                    if ((aa != bb && bb != dd) && (aa != dd && bb != cc))
                      r = ddfVal(cc) / 4;
                    else
                      r = 0;
 
                  } else if ((cc != 1 && dd != 1) && (aa == 1 || bb == 1))
                    r = dfVal(cc) * aF(aa, bb) / 2;
                  else if ((cc == 1 || dd == 1) && (aa != 1 && bb != 1)) {
 
                    if ((cc == 2 && dd == 1) || (cc == 1 && dd == 2))
                      r = (
 
                              dfVal(cc)
 
                              - (fVal(cc) - fVal(dd)) * aF(cc, dd)
 
                                  ) *
                          aF(cc, dd);
 
                    else
                      r = 0;
 
                  } else
                    r = 0;
 
                  if (r)
                    d2MType1[dd][get_sym_index(aa, cc, Dim)](i, j, k, l) +=
                        r * aS[get_sym_index(aa, bb, Dim)](i, j, k, l);
                }
            }
        }
      }
 
    if constexpr (NB == 1)
      for (auto aa = 0; aa != Dim; ++aa) {
        for (auto bb = 0; bb != Dim; ++bb) {
          if (aa != bb) {
            for (auto cc = 0; cc != Dim; ++cc) {
              for (auto dd = 0; dd != Dim; ++dd) {
                if (cc != dd) {
                  if ((bb != dd) && (aa != dd && bb != cc)) {
                    const double r = ddfVal(cc) / 4;
                    d2MType1[dd][get_sym_index(aa, cc, Dim)](i, j, k, l) +=
                        r * aS[get_sym_index(aa, bb, Dim)](i, j, k, l);
                  }
                }
              }
            }
          }
        }
      }
 
    using THIS = EigenMatrixImp<T1, T2, NB, Dim>;
    using T3 = FTensor::Ddg<V, Dim, Dim>;
 
    T3 t_diff_A;
    GetDiffDiffMatImpl<THIS, V, T3, T>(*this, t_diff_A, t_S)
        .set(Number<Dim>(), Number<Dim>(), Number<Dim>(), Number<Dim>());
    return t_diff_A;
  }

◆ getDiffMat()

template<typename T1 , typename T2 , int NB, int Dim>

auto EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::getDiffMat	(	Fun	f,
		Fun	d_f
	)

inline

Get derivative of matrix.

\[ P_{ijkl} = \frac{\partial B_{ij}}{\partial A_{kl}} \]

Parameters

t_val	eiegn values vector
t_vec	eiegn vectors matrix
f	function
d_f	directive of function

Returns: auto derivatives, forth order tensor with minor simetries

Definition at line 1032 of file MatrixFunctionTemplate.hpp.

                                         {
 
    for (auto aa = 0; aa != Dim; ++aa)
      fVal(aa) = f(tVal(aa));
 
    for (auto aa = 0; aa != Dim; ++aa)
      dfVal(aa) = d_f(tVal(aa));
 
    for (auto aa = 0; aa != Dim; ++aa)
      for (auto bb = 0; bb != aa; ++bb) {
        aF(aa, bb) = 1 / (tVal(aa) - tVal(bb));
        aF(bb, aa) = -aF(aa, bb);
        aF2(aa, bb) = aF(aa, bb) * aF(aa, bb);
      }
 
    using T3 = FTensor::Ddg<V, Dim, Dim>;
    T3 t_diff_A;
    GetDiffMatImpl<EigenMatrixImp<T1, T2, NB, Dim>, V, T3>(*this, t_diff_A)
        .set(NumberDim(), NumberDim(), NumberDim(), NumberDim());
    return t_diff_A;
  }

◆ getMat()

template<typename T1 , typename T2 , int NB, int Dim>

auto EigenMatrix::EigenMatrixImp< T1, T2, NB, Dim >::getMat ( Fun f )

inline

Get matrix.

\[ \mathbf{B} = f(\mathbf{A}) \]

\[ B_{ij} = \sum_{a}^3 f(\lambda^a) n^a_i n^a_j \]

where \(a\) is eigen value number.

Parameters

t_val	eiegn values vector
t_vec	eigen vectors matrix
f	function

Returns: auto function symmetric tensor rank two

Definition at line 1006 of file MatrixFunctionTemplate.hpp.

                            {
 
    for (auto aa = 0; aa != Dim; ++aa)
      fVal(aa) = f(tVal(aa));
 
    using T3 =
        FTensor::Tensor2_symmetric<typename std::remove_const<V>::type, Dim>;
    T3 t_A;
    GetMatImpl<EigenMatrixImp<T1, T2, NB, Dim>, V, T3>(*this, t_A)
        .set(NumberDim(), NumberDim());
    return t_A;
  }