#ifndef DUNE_ELASTICITY_HELPERS_HH
#define DUNE_ELASTICITY_HELPERS_HH

#include <cmath>

#include <dune/common/fmatrix.hh>

namespace Dune {

    namespace Elasticity {

        /** \brief Compute nonlinear strain tensor from a displacement gradient.
         *
         *  \param grad The gradient of the direction in which the linearisation is computed.
         *  \param strain The tensor to store the strain in. 
         */
        template <int dim>
        void computeNonlinearStrain(const Dune::FieldMatrix<double, dim, dim>& grad, SymmetricTensor<dim>& strain) {
            strain = 0;
            for (int i=0; i<dim ; ++i) {
                strain(i,i) +=grad[i][i];

                for (int k=0;k<dim;++k)
                    strain(i,i) += 0.5*grad[k][i]*grad[k][i];

                for (int j=i+1; j<dim; ++j) {
                    strain(i,j) += 0.5*(grad[i][j] + grad[j][i]);
                    for (int k=0;k<dim;++k)
                        strain(i,j) += 0.5*grad[k][i]*grad[k][j];
                }
            }
        }

        /** \brief Compute linearised strain at the identity from a given displacement gradient. 
         *
         *  \param grad The gradient of the direction in which the linearisation is computed.
         *  \param strain The tensor to store the strain in. 
         */
        template <int dim>
        void computeLinearisedStrain(const Dune::FieldMatrix<double, dim, dim>& grad, SymmetricTensor<dim>& strain) {
            for (int i=0; i<dim ; ++i)
            {
                strain(i,i) = grad[i][i];
                for (int j=i+1; j<dim; ++j)
                    strain(i,j) = 0.5*(grad[i][j] + grad[j][i]);
            }
        }


        /** \brief Compute linearised strain at a prescribed configuration in the direction of a given displacement gradient. 
         *
         *  \param grad The gradient of the direction in which the linearisation is computed.
         *  \param conf The deformation gradient(!) of the configuration at which the linearisation is evaluated.
         *  \param strain The tensor to store the strain in. 
         */
        template <int dim>
        void computeLinearisedStrain(const Dune::FieldMatrix<double, dim, dim>& grad, const Dune::FieldMatrix<double, dim, dim>& conf,
                                     SymmetricTensor<dim>& strain) {
            strain = 0;
            for (int i=0;i<dim; i++)
                for (int j=i;j<dim;j++)
                    for (int k=0;k<dim;k++)
                        strain(i,j)+= grad[k][i]*conf[k][j]+conf[k][i]*grad[k][j];
            strain *= 0.5;
        }

        /** \brief The penalty term \f$ \Gamma \f$ */
        double Gamma(double x) {
            return x*x - std::log(x);
        }
        
        /** \brief First derivative of the penalty term \f$ \Gamma \f$ */
        double Gamma_x(double x) {
            return 2*x - 1/x;
        }
        
        /** \brief Second derivative of the penalty term \f$ \Gamma \f$ */
        double Gamma_xx(double x) {
            return 2 + 1/(x*x);
        }
        
        double delta(int i, int j) { return i==j? 1.0: 0.0; }
        
        template <int dim>
        double E_du(const Dune::FieldMatrix<double,dim,dim>& u, int i, int j, int p, int q) 
        {
            double result = delta(i,p)*delta(j,q) + delta(j,p)*delta(i,q);
            
            result += ( delta(i,q)*u[p][j] + u[p][i]*delta(j,q) );
            return 0.5*result;
        }
        
        // template <int i, int j, int p, int q, int r, int s>
        // static inline double e_dudu() 
        // {
        //     return 0.5 * Delta<p,r>::delta
        //         * ( Delta<i,q>::delta*Delta<j,s>::delta + Delta<i,s>::delta*Delta<j,q>::delta );
        // }
        
        static inline double e_dudu(int i, int j, int p, int q, int r, int s) 
        {
            return 0.5 * delta(p,r)* ( delta(i,q)*delta(j,s) + delta(i,s)*delta(j,q) );
        }
        
        // The Kronecker delta function
        template <int i, int j>
        struct Delta {
            static int const delta = 0;
        };
        
        template <int i>
        struct Delta<i,i> {
            static int const delta = 1;
        };
        
        // \parder {det I + u}{ u_pq }
        template <int p, int qq>
        double det_du_tmpl(const Dune::FieldMatrix<double,3,3>& u) {
            int const q = qq+3;
            
            return (Delta<(p+1)%3,(q+1)%3>::delta + u[(p+1)%3][(q+1)%3])
                * (Delta<(p+2)%3,(q+2)%3>::delta+u[(p+2)%3][(q+2)%3]) 
                - (Delta<(p+1)%3,(q-1)%3>::delta+u[(p+1)%3][(q-1)%3])
                * (Delta<(p+2)%3,(q-2)%3>::delta+u[(p+2)%3][(q-2)%3]);
        }
         
        /** \brief Compute linearization of the determinant of a deformation gradient
         *
         *  \param u The displacement gradient(!) at which the determinant is evaluated
         */
        void linearizedDefDet(const Dune::FieldMatrix<double,3,3>& u, Dune::FieldMatrix<double,3,3>& linDet) {

            linDet = 0;
            for (int i=0; i<2; i++)                       
                for (int j=(i+1)%3; j<3; j++) {
                    int k=(j+1)%3 + (j+1)/3;          
                    linDet[i][j] = u[j][k]*u[k][i] - u[j][i]*(1+u[k][k]);
                    linDet[j][i] = u[k][j]*u[i][k] - u[i][j]*(1+u[k][k]);
            }    

            // the diagonal parts
            for (int i=0; i<3; i++)
                linDet[i][i] = (u[(i+1)%3][(i+1)%3]+1)*(u[(i+2)%3][(i+2)%3]+1)-u[(i+1)%3][(i+2)%3]*u[(i+2)%3][(i+1)%3]
        }
        
        /** \brief Compute linearization of the determinant of a deformation gradient
         *
         *  \param u The displacement gradient(!) at which the determinant is evaluated
         */
        void linearizedDefDet(const Dune::FieldMatrix<double,2,2>& u, Dune::FieldMatrix<double,2,2>& linDet) {

            linDet = 0;
            for (int i=0; i<2; i++) {
                linDet[i][i] = 1+u[(i+1)%2][(i+1)%2];
                linDet[i][(i+1)%2] = -u[(i+1)%2][i];
            }
        }

        // \parder {det I + u}{ u_pq }
        template <int p, int q>
        double det_du_tmpl(const Dune::FieldMatrix<double,2,2>& u) {
            return Delta<p,q>::delta * (1+u[(p+1)%2][(q+1)%2]) 
                - Delta<p,(q+1)%2>::delta * u[(p+1)%2][(q+1)%2];
            
            // p = 0, q = 0;  (1+u[1][1]) 
            // p = 0, q = 1; -u[1][0] 
            // p = 1, q = 0; -u[0][1]
            // p = 1, q = 1;  (1+u[0][0])
            //
        }
        
        // template <int p, int qq, int r, int s>
        // double det_dudu_tmpl(const Dune::FieldMatrix<double,3,3>& u) {
        //     int const q = qq+3;
        
        //     return
        //         Delta<(p+1)%3,r>::delta*Delta<(q+1)%3,s>::delta * (Delta<(p+2)%3,(q+2)%3>::delta+u[(p+2)%3][(q+2)%3])
        //         + Delta<(p+2)%3,r>::delta*Delta<(q+2)%3,s>::delta * (Delta<(p+1)%3,(q+1)%3>::delta+u[(p+1)%3][(q+1)%3])
        //         - Delta<(p+1)%3,r>::delta*Delta<(q-1)%3,s>::delta * (Delta<(p+2)%3,(q-2)%3>::delta+u[(p+2)%3][(q-2)%3])
        //         - Delta<(p+2)%3,r>::delta*Delta<(q-2)%3,s>::delta * (Delta<(p+1)%3,(q-1)%3>::delta+u[(p+1)%3][(q-1)%3]);
        // }
        
        double det_dudu_tmpl(const Dune::FieldMatrix<double,3,3>& u, int p, int qq, int r, int s) {
            int const q = qq+3;
            
            return
                delta((p+1)%3,r) * delta((q+1)%3,s) * (delta((p+2)%3,(q+2)%3)+u[(p+2)%3][(q+2)%3])
                + delta((p+2)%3,r)*delta((q+2)%3,s) * (delta((p+1)%3,(q+1)%3)+u[(p+1)%3][(q+1)%3])
                - delta((p+1)%3,r)*delta((q-1)%3,s) * (delta((p+2)%3,(q-2)%3)+u[(p+2)%3][(q-2)%3])
                - delta((p+2)%3,r)*delta((q-2)%3,s) * (delta((p+1)%3,(q-1)%3)+u[(p+1)%3][(q-1)%3]);
        }
        
        // parder det (I + u) partial u_pq \partial u_rs
        //template <int p, int q, int r, int s>
        double det_dudu_tmpl(const Dune::FieldMatrix<double,2,2>& u, int p, int q, int r, int s) {
            return delta(p,q)*delta(r,(p+1)%2)*delta(s,(q+1)%2)
                - delta(p,(q+1)%2)*delta(r,(p+1)%2)*delta(s,(q+1)%2);
            
            // 0 0 0 0   -->  0  - 0
            // 0 0 0 1   -->  0  - 0
            // 0 0 1 0   -->  0  - 0
            // 0 0 1 1   -->  1  - 0  = 1
            // 0 1 0 0   -->  0  - 0
            // 0 1 0 1   -->  0  - 0
            // 0 1 1 0   -->  0  - 1  = -1
            // 0 1 1 1   -->  0  - 0
            // 1 0 0 0   -->  0  - 0
            // 1 0 0 1   -->  0  - 1  = -1
            // 1 0 1 0   -->  0  - 0
            // 1 0 1 1   -->  0  - 0
            // 1 1 0 0   -->  1  - 0  = 1
            // 1 1 0 1   -->  0  - 0
            // 1 1 1 0   -->  0  - 0
            // 1 1 1 1   -->  0  - 0
        }
        
        double det_val(const Dune::FieldMatrix<double,3,3>& u) {
            return (1+u[0][0])*(1+u[1][1])*(1+u[2][2])
                + u[0][1]*u[1][2]*u[2][0]
                + u[0][2]*u[1][0]*u[2][1]
                - u[2][0]*(1+u[1][1])*u[0][2]
                - u[2][1]*u[1][2]*(1+u[0][0])
                - (1+u[2][2])*u[1][0]*u[0][1];
        }
        
        double det_val(const Dune::FieldMatrix<double,2,2>& u) {
            return (1+u[0][0])*(1+u[1][1]) - u[1][0]*u[0][1];
        }
        
        template <int dim>
        double E_val(const Dune::FieldMatrix<double,dim,dim>& u, int i, int j) 
        {
            double result = u[i][j] + u[j][i];
            
            for (int k=0; k<dim; ++k)
                result += u[k][i]*u[k][j];  
            return 0.5*result;
        }
        
        

    }

}

#endif