/****************************************************************************
*
*                                 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 <iostream>
#include <string>
#include <stdlib.h>
#include "newtonian.h"


using namespace std;
using namespace muesli;


newtonianMaterial::newtonianMaterial(const std::string& name)
:
fluidMaterial(name),
mu(0.0), lambda(0.0), _incompressible(true), theEOS(_pg)
{}




newtonianMaterial::newtonianMaterial(const std::string& name,
                                       const materialProperties& cl) :
  fluidMaterial(name, cl),
  mu(0.0), lambda(0.0), _incompressible(true), theEOS(_pg)
{

    muesli::assignValue(cl, "mu"    , mu);
    muesli::assignValue(cl, "lambda", lambda);

    string eosname;
    if (muesli::assignValue(cl, "eos", eosname))
    {
        _incompressible = false;
        if (eosname.compare("_pg"))
        {   
            theEOS = _pg;
            double gamma = 1.4;
            muesli::assignValue(cl, "gamma", gamma);
            eosConst.push_back(gamma);
        }
        else
        {   
            std::cout << "Equation of state not found" << std::endl;
            exit(-1);
        }
    }
}




bool newtonianMaterial::check() const
{
    bool ret = true;
    if (mu <= 0.0)
    {
        std::cout << "Error in newtonianMaterial material. Zero viscosity.";
        ret = false;
    }

    return ret;
}




muesli::fluidMP* newtonianMaterial::createMaterialPoint() const
{
    muesli::fluidMP* mp = new newtonianMP(*this);
    return mp;
}




double newtonianMaterial::getProperty(const propertyName p) const
{
    double ret=0.0;
    double rho = density();

    // scan all the possible data
    switch(p)
    {
        case PR_MU:         ret = mu;           break;
        case PR_LAMBDA:     ret = lambda;       break;
        case PR_NU:         ret = mu/rho;       break;
        case PR_CS:         ret = sqrt(mu/rho); break;

        default:
            std::cout << "Error in newtonianMaterial. Property not defined";
    }
    return ret;

}




double newtonianMaterial::kinematicViscosity() const
{
    return mu/density();
}




/*! this function is always called once the material is defined, so apart from
 printing its information, we take the opportunity to clean up some of its
 data, in particular, setting all the possible constants
 */
void newtonianMaterial::print(std::ostream &of) const
{
    if (_incompressible)
        of  << "\n   Incompressible Newtonian fluid";
    else
        of  << "\n   Compressible Newtonian fluid";

    of  << "\n   Kinematic viscosity  nu     : " << mu/density()
        << "\n   Dynamic viscosity    mu     : " << mu
        << "\n   2nd viscosity        lambda : " << lambda
        << "\n   Density:                    : " << density();

    if (density() > 0.0)
        of  << "\n   Wave velocity   c_s         : " << sqrt(mu/density());

    of  << "\n";
}




void newtonianMaterial::setRandom()
{
    setDensity(muesli::randomUniform(1.0, 100.0));
    mu      = muesli::randomUniform(0.0, 1.0);
    lambda  = muesli::randomUniform(0.0, 1.0);
    _incompressible = true;
}




double newtonianMaterial::waveVelocity() const
{
    return sqrt(mu/density());
}




bool newtonianMaterial::test(std::ostream &of)
{
    bool isok = true;
    setRandom();

    fluidMP* p = this->createMaterialPoint();
    isok = p->testImplementation(of);
    delete p;

    return isok;
}




newtonianMP::newtonianMP(const newtonianMaterial &m) :
    fluidMP(const_cast<newtonianMaterial &>(m)),
    theNewtonianMaterial(&m),
    d_c(m.density())
{
}




void newtonianMP::CauchyStress(istensor &sigma) const
{
    const double mu     = theNewtonianMaterial->mu;
    const double lambda = theNewtonianMaterial->lambda;

    istensor D = istensor::symmetricPartOf(gradu_c);
    sigma = 2.0 * mu * D;

    sigma += (-p_c + ( lambda - 2.0/3.0 * mu ) * gradu_c.trace()) * istensor::identity();
}




void newtonianMP::commitCurrentState()
{
    t_n     = t_c;
    gradu_n = gradu_c;
    p_n     = p_c;
    d_n     = d_c;
}




void newtonianMP::contractWithDeviatoricTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    const double&  mu     = theNewtonianMaterial->mu;
    const double&  lambda = theNewtonianMaterial->lambda;

    T = mu * v1.dot(v2)         * itensor::identity()
    +   mu                      * itensor::dyadic(v2, v1)
    + (lambda -  2.0/3.0*mu )   * itensor::dyadic(v1, v2);
}




double newtonianMP::density() const
{
    return d_c;
}




double newtonianMP::deviatoricEnergy() const
{
    istensor D;
    deviatoricStress(D);
    return D.contract(gradu_c);
}




void newtonianMP::deviatoricStress(istensor &sigma) const
{
    const double mu = theNewtonianMaterial->mu;

    istensor D = istensor::symmetricPartOf(gradu_c);
    sigma = 2.0 * mu * D;
    sigma -= 2.0/3.0 * mu  * gradu_c.trace() * istensor::identity();
}




double newtonianMP::dissipatedEnergy() const
{
    return 0.0;
}




double newtonianMP::pressure() const
{
    return p_c;
}




double newtonianMP::storedEnergy() const
{
    return volumetricEnergy() + deviatoricEnergy();
}




void newtonianMP::tangentTensor(itensor4& C) const
{
    C.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++)
                {
                    if ((i==k) && (j==l)) C(i,j,k,l) += theNewtonianMaterial->mu;
                    if ((i==l) && (j==k)) C(i,j,k,l) += theNewtonianMaterial->mu;
                    if ((i==j) && (k==l)) C(i,j,k,l) += theNewtonianMaterial->lambda - (2.0/3.0)*theNewtonianMaterial->mu;
                }
}




void newtonianMP::updateCurrentState(const double theTime, const itensor& gradu)
{
    t_c     = theTime;
    gradu_c = gradu;
    d_c     = theNewtonianMaterial->density();
}





void newtonianMP::updateCurrentState(const double theTime, const itensor& gradu, const double pressure, const double rho)
{
    t_c     = theTime;
    gradu_c = gradu;
    p_c     = pressure;
    d_c     = rho;
}




void newtonianMP::updateCurrentState(const double theTime, const itensor& gradu, const double pressure)
{
    t_c     = theTime;
    gradu_c = gradu;
    p_c     = pressure;
    d_c     = theNewtonianMaterial->density();
}




double newtonianMP::volumetricEnergy() const
{
    istensor V;
    volumetricStress(V);
    return V.contract(gradu_c);
}




void newtonianMP::volumetricStress(istensor &sigma) const
{
    const double lambda = theNewtonianMaterial->lambda;
    sigma = (-p_c + lambda * gradu_c.trace()) * istensor::identity();
}

