principalInvariants.hh

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                             */
/*   This file is part of the C++-library RFFGen.                              */
/*   Copyright 2015 Lars Lubkoll                                               */
/*                                                                             */
/*   RFFGen is free software: you can redistribute it and/or modify            */
/*   it under the terms of the GNU General Public License as published by      */
/*   the Free Software Foundation, either version 3 of the License, or         */
/*   (at your option) any later version.                                       */
/*                                                                             */
/*   RFFGen is distributed in the hope that it will be useful,                 */
/*   but WITHOUT ANY WARRANTY; without even the implied warranty of            */
/*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             */
/*   GNU General Public License for more details.                              */
/*                                                                             */
/*   You should have received a copy of the GNU General Public License         */
/*   along with RFFGen.  If not, see <http://www.gnu.org/licenses/>.           */
/*                                                                             */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#ifndef RFFGEN_LINEAR_ALGEBRA_PRINCIPAL_INVARIANTS_HH
#define RFFGEN_LINEAR_ALGEBRA_PRINCIPAL_INVARIANTS_HH

#include <type_traits>

#include "dimension.hh"
#include "trace.hh"
#include "cofactor.hh"
#include "determinant.hh"
#include "shiftedInvariant.hh"
#include "../Util/base.hh"

namespace RFFGen
{
  namespace Concepts { template <class> struct SymmetricMatrixConceptCheck; }
  namespace LinearAlgebra
  {
    namespace Detail
    {
      template <class Matrix>
      auto computeCofactors(const Matrix& A, std::integral_constant<int,2>)
      {
        return computeCofactor<0,0>(A) + computeCofactor<1,1>(A);
      }

      template <class Matrix>
      auto computeCofactors(const Matrix& A, std::integral_constant<int,3>)
      {
        return computeCofactor<0,0>(A) + computeCofactor<1,1>(A) + computeCofactor<2,2>(A);
      }

      template < class Matrix >
      auto computeCofactors(const Matrix& A, std::integral_constant<int,-1>)
      {
        if(rows(A) == 2) return computeCofactors(A,std::integral_constant<int,2>());
        /*if(rows(A) == 3)*/ return computeCofactors(A,std::integral_constant<int,3>());
      }

      template <int row, class Matrix>
      auto symmetricCofactorDerivative(const Matrix& A, const Matrix& dA)
      {
        return computeCofactorDirectionalDerivative<row,row>(A,dA) + computeCofactorDirectionalDerivative<row,row>(dA,A);
      }


      template <class Matrix>
      auto computeSymmetricCofactorDerivatives(const Matrix& A, const Matrix& B, std::integral_constant<int,2>)
      {
        return symmetricCofactorDerivative<0>(A,B) + symmetricCofactorDerivative<1>(A,B);
      }

      template <class Matrix>
      auto computeSymmetricCofactorDerivatives(const Matrix& A, const Matrix& B, std::integral_constant<int,3>)
      {
        return symmetricCofactorDerivative<0>(A,B) + symmetricCofactorDerivative<1>(A,B) + symmetricCofactorDerivative<2>(A,B);
      }


      template < class Matrix >
      auto computeSymmetricCofactorDerivatives(const Matrix& A, const Matrix& B, std::integral_constant<int,-1>)
      {
        if(rows(A) == 2) return computeSymmetricCofactorDerivatives(A,B,std::integral_constant<int,2>());
        /*if(rows(A) == 3)*/ return computeSymmetricCofactorDerivatives(A,B,std::integral_constant<int,3>());

      }
    }
    template <class Matrix>
    using FirstPrincipalInvariant = Trace<Matrix>;

    template <class Matrix, class = Concepts::SymmetricMatrixConceptCheck<Matrix> >
    class SecondPrincipalInvariant : Base
    {
      using DimIdentifier = std::integral_constant<int,dimension<Matrix>()>;
    public:
      SecondPrincipalInvariant() = default;

      SecondPrincipalInvariant(const Matrix& A) { update(A); }

      void update(const Matrix& A)
      {
        A_ = A;
        resultOfD0 = Detail::computeCofactors(A,dimIdentifier);
      }

      auto d0() const
      {
        return resultOfD0;
      }

      template <int,int>
      auto d1(const Matrix& dA1) const
      {
        return Detail::computeSymmetricCofactorDerivatives(A_,dA1,dimIdentifier);
      }

      template <int,int>
      auto d2(const Matrix& dA1, const Matrix& dA2) const
      {
        return Detail::computeSymmetricCofactorDerivatives(dA1,dA2,dimIdentifier);
      }

    private:
      Matrix A_;
      std::decay_t< at_t<Matrix> > resultOfD0 = 0.;
      DimIdentifier dimIdentifier;
    };

    template <class Matrix>
    using ThirdPrincipalInvariant = Determinant<Matrix>;


    template < class Matrix , int offset = dimension<Matrix>() >
    using ShiftedFirstPrincipalInvariant  = ShiftedInvariant< FirstPrincipalInvariant<Matrix> , offset >;

    template < class Matrix , int offset = dimension<Matrix>() >
    using ShiftedSecondPrincipalInvariant = ShiftedInvariant< SecondPrincipalInvariant<Matrix> , offset >;

    template < class Matrix >
    using ShiftedThirdPrincipalInvariant  = ShiftedInvariant< ThirdPrincipalInvariant<Matrix>   , 1 >;
  }
}

#endif // RFFGEN_LINEAR_ALGEBRA_PRINCIPAL_INVARIANTS_HH