/****************************************************************************
*
*                                 M U E S L I   v 1.8
*
*
*     Copyright 2020 IMDEA Materials Institute, Getafe, Madrid, Spain
*     Contact: muesli.materials@imdea.org
*     Author: Ignacio Romero (ignacio.romero@imdea.org)
*
*     This file is part of MUESLI.
*
*     MUESLI 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.
*
*     MUESLI 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 MUESLI.  If not, see <http://www.gnu.org/licenses/>.
*
****************************************************************************/



#include "fluid.h"
#include <iostream>

using namespace muesli;


fluidMaterial::fluidMaterial(const std::string& name) :
material(name)
{
}




fluidMaterial::fluidMaterial(const std::string& name,
                               const materialProperties& cl)  :
  material(name, cl)
{
}



fluidMP::fluidMP(const fluidMaterial &m) :
theFluidMaterial(m)
{
}



materialState fluidMP::getConvergedState() const
{
    materialState state;

    return state;
}




materialState fluidMP::getCurrentState() const
{
    materialState state;

    return state;
}




double fluidMP::pressure(eos theEOS, std::vector<double> eosConst, double rho, double e) const
{
    double p = 0.0;
    if (theEOS == _pg)
    {       
        double gamma = eosConst[0];     
        p = rho*e*(gamma-1.0);
    }
    return p;
}




double fluidMP::energy(eos theEOS, std::vector<double> eosConst, double rho, double p) const
{
    double e = 0.0;
    if (theEOS == _pg)
    {       
        double gamma = eosConst[0];
        e = p/(rho*(gamma-1.0));
    }
    return e;
}




void fluidMP::contractWithTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    itensor4 C;
    tangentTensor(C);

    T.setZero();
    for (unsigned i=0; i<3; i++)
        for (unsigned j=0; j<3; j++)
            for (unsigned k=0; k<3; k++)
                for (unsigned l=0; l<3; l++)
                    T(i,k) += C(i,j,k,l)*v1[j]*v2[l];
}




/* let Pdev the 4th order tensor that projects sym tensors to deviatoric sym tensors
   Compute Cdev = Pdev^T C Pdev, then
   compute the contraction (Cdev)_{ijkl} v_j w_l
 */
void fluidMP::contractWithDeviatoricTangent(const ivector &v, const ivector &w, itensor &T) const
{
    itensor4 C;
    tangentTensor(C);

    double Cpp[3][3];
    for (unsigned i=0; i<3; i++)
        for (unsigned j=0; j<3; j++)
        {
            Cpp[i][j] = 0.0;
            for (unsigned p=0; p<3; p++)
                Cpp[i][j] += C(p,p,i,j);
        }

    double Cppqq=0.0;
    for (unsigned q=0; q<3; q++)
    {
        Cppqq += Cpp[q][q];
    }


    T.setZero();
    const double th = 1.0/3.0;
    for (unsigned i=0; i<3; i++)
    {
        for (unsigned k=0; k<3; k++)
        {
            for (unsigned j=0; j<3; j++)
            {
                for (unsigned l=0; l<3; l++)
                {
                    T(i,k) += C(i,j,k,l)*v[j]*w[l];
                }

                T(i,k) -= th*Cpp[k][j]*v[i]*w[j]
                        + th*Cpp[i][j]*v[j]*w[k];
            }
            T(i,k) += th*th*Cppqq *v[i]*w[k];
        }
    }
}




bool fluidMP::testImplementation(std::ostream& os) const
{

    bool isok = true;
    // set a random update in the material
    istensor eps;
    eps.setZero();
    double p = 0.0;
    const_cast<fluidMP*>(this)->updateCurrentState(0.0, eps, p);
    const_cast<fluidMP*>(this)->commitCurrentState();

    double tn1 = muesli::randomUniform(0.1,1.0);
    eps.setRandom();
    double pn1 = muesli::randomUniform(0.1,1.0);
    const_cast<fluidMP*>(this)->updateCurrentState(tn1, eps, pn1);

    // tangent
    itensor4 tg;
    tangentTensor(tg);

    // compare tensor c with derivative of stress
    if ((true))
    {
        // numeric C
        itensor4 nC;
        nC.setZero();

        // numerical differentiation sigma
        istensor dsigma, sigmap1, sigmap2, sigmam1, sigmam2;
        double   inc = 1.0e-3;

        for (unsigned i=0; i<3; i++)
        {
            for (unsigned j=i; j<3; j++)
            {
                double  original = eps(i,j);

                eps(i,j) = eps(j,i) = original + inc;
                const_cast<fluidMP*>(this)->updateCurrentState(tn1, eps, pn1);
                CauchyStress(sigmap1);

                eps(i,j) = eps(j,i) = original + 2.0*inc;
                const_cast<fluidMP*>(this)->updateCurrentState(tn1, eps, pn1);
                CauchyStress(sigmap2);

                eps(i,j) = eps(j,i) = original - inc;
                const_cast<fluidMP*>(this)->updateCurrentState(tn1, eps, pn1);
                CauchyStress(sigmam1);

                eps(i,j) = eps(j,i) = original - 2.0*inc;
                const_cast<fluidMP*>(this)->updateCurrentState(tn1, eps, pn1);
                CauchyStress(sigmam2);

                // fourth order approximation of the derivative
                dsigma = (-sigmap2 + 8.0*sigmap1 - 8.0*sigmam1 + sigmam2)/(12.0*inc);

                if (i != j) dsigma *= 0.5;


                for (unsigned k=0; k<3; k++)
                    for (unsigned l=0; l<3; l++)
                    {
                        nC(k,l,i,j) = dsigma(k,l);
                        nC(k,l,j,i) = dsigma(k,l);
                    }

                eps(i,j) = original;
                eps(j,i) = original;
                const_cast<fluidMP*>(this)->updateCurrentState(tn1, eps, pn1);
            }
        }

        // relative error less than 0.01%
        double error = 0.0;
        double norm = 0.0;
        for (unsigned i=0; i<3; i++)
            for (unsigned j=0; j<3; j++)
                for (unsigned k=0; k<3; k++)
                    for (unsigned l=0; l<3; l++)
                    {
                        error += pow(nC(i,j,k,l)-tg(i,j,k,l),2);
                        norm  += pow(tg(i,j,k,l),2);
                    }
        error = sqrt(error);
        norm = sqrt(norm);
        isok = (error/norm < 1e-4);

        os << "\n   1. Comparing tensor C with DStress.";
        if (isok)
        {
            os << " Test passed.";
        }
        else
        {
            os << "\n      Test failed.";
            os << "\n      Relative error in DWeff computation: " <<  error/norm;
        }

    }

    return isok;
}
