/****************************************************************************
*
*                                 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 "conductor.h"


using namespace std;
using namespace muesli;



conductorMaterial::conductorMaterial(const std::string& name)
:
material(name)
{

}




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

}




bool conductorMaterial::check() const
{
    bool ret = true;
    return ret;
}




void conductorMaterial::setRandom()
{
    setDensity(muesli::randomUniform(1.0, 10.0));
}




conductorMP::conductorMP(const conductorMaterial& mat) :
theConductor(&mat),
time_n(0.0),
time_c(0.0)
{
    gradT_n.setZero();
    gradT_c.setZero();
}




void conductorMP::commitCurrentState()
{
    time_n  = time_c;
    temp_n  = temp_c;
    gradT_n = gradT_c;
}




void conductorMP::resetCurrentState()
{
    time_c  = time_n;
    temp_c  = temp_n;
    gradT_c = gradT_n;
}




void conductorMP::setRandom()
{
    gradT_n.setRandom();
    gradT_c = gradT_n;
}




bool conductorMP::testImplementation(std::ostream& os) const
{
    bool isok = true;

    double tn1 = muesli::randomUniform(0.1,1.0);

    ivector gradT; gradT.setRandom();
    double temp = muesli::randomUniform(280, 350);
    const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);

    // programmed heatflux and conductivity
    ivector q;  heatflux(q);
    istensor K; conductivity(K);
    ivector qp; heatfluxDerivative(qp);

    // compare heat flux with (minus) derivative of energy
    if (true)
    {
        // numerical differentiation heatflux
        ivector numFlux;
        numFlux.setZero();
        const double inc = 1.0e-3;

        for (size_t i=0; i<3; i++)
        {
            double original = gradT(i);

            gradT(i) = original + inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            double Wp1 = thermalEnergy();

            gradT(i) = original + 2.0*inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            double Wp2 = thermalEnergy();

            gradT(i) = original - inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            double Wm1 = thermalEnergy();

            gradT(i) = original - 2.0*inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            double Wm2 = thermalEnergy();


            // fourth order approximation of the derivative
            double der = - (-Wp2 + 8.0*Wp1 - 8.0*Wm1 + Wm2)/(12.0*inc);
            numFlux(i) = der;

            gradT(i) = original;
        }

        // relative error less than 0.01%
        ivector error = numFlux - q;
        isok = (error.norm()/q.norm() < 1e-4);

        os << "\n   1. Comparing heat flux with derivative of thermal energy.";

        if (!isok)
        {
            os << "\n Relative error in q computation %e. Test failed." <<  error.norm()/q.norm();
            os << "\n   Heat flux:\n" << q;
            os << "\n   Numeric heat flux:\n" << numFlux;
        }
        else
            os << " Test passed.";
    }
    

    // compare conductivity tensor with (minus) derivative of heat flux wrt temp gradient
    if (true)
    {
        istensor numK;
        numK.setZero();
        const double inc = 1.0e-3;

        for (size_t i=0; i<3; i++)
        {
            double original = gradT(i);

            gradT(i) = original + inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            ivector Qp1; heatflux(Qp1);

            gradT(i) = original + 2.0*inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            ivector Qp2; heatflux(Qp2);

            gradT(i) = original - inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            ivector Qm1; heatflux(Qm1);

            gradT(i) = original - 2.0*inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            ivector Qm2; heatflux(Qm2);

            // fourth order approximation of the derivative
            ivector der = - (-Qp2 + 8.0*Qp1 - 8.0*Qm1 + Qm2)/(12.0*inc);
            numK(0,i) = der(0);
            numK(1,i) = der(1);
            numK(2,i) = der(2);

            gradT(i) = original;
        }

        // relative error less than 0.01%
        istensor error = numK - K;
        isok = (error.norm()/q.norm() < 1e-4);

        os << "\n   2. Comparing conductivity tensor with derivative of heat flux.";

        if (!isok)
        {
            os << "\n Relative error in q computation %e. Test failed." <<  error.norm()/q.norm();
            os << "\n   Conductivity flux:\n" << K;
            os << "\n   Numeric conductivity :\n" << numK;
        }
        else
            os << " Test passed.";
    }

    // compare qprime with derivative of heat flux wrt temp
    if (true)
    {
        ivector numQ;
        numQ.setZero();
        const double inc = 1.0e-3;

        {
            double original = temp;

            temp = original + inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            ivector Qp1; heatflux(Qp1);

            temp = original + 2.0*inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            ivector Qp2; heatflux(Qp2);

            temp = original - inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            ivector Qm1; heatflux(Qm1);

            temp = original - 2.0*inc;
            const_cast<conductorMP*>(this)->updateCurrentState(tn1, temp, gradT);
            ivector Qm2; heatflux(Qm2);

            // fourth order approximation of the derivative
            ivector der = (-Qp2 + 8.0*Qp1 - 8.0*Qm1 + Qm2)/(12.0*inc);
            numQ = der;

            temp = original;
        }

        // relative error less than 0.01%
        ivector error = numQ - qp;
        isok = (error.norm()/qp.norm() < 1e-4);

        os << "\n   3. Comparing qprime with derivative of heat flux.";

        if (!isok)
        {
            os << "\n Relative error in q computation %e. Test failed." <<  error.norm()/q.norm();
            os << "\n   Q prime: " << qp;
            os << "\n   Numeric derivative : " << numQ;
        }
        else
            os << " Test passed.";
    }


    return isok;
}




void conductorMP::updateCurrentState(double theTime, double temp, const ivector& gradT)
{
    time_c  = theTime;
    temp_c  = temp;
    gradT_c = gradT;
}




fourierMaterial::fourierMaterial(const std::string& name,
                                   const materialProperties& cl)
:
  conductorMaterial(name, cl),
  _conductivity(0.0),
  _capacity(0.0)
{
    muesli::assignValue(cl, "conductivity", _conductivity);
    muesli::assignValue(cl, "capacity", _capacity);
}




fourierMaterial::fourierMaterial(const std::string& name, const double k, const double rho)
:
conductorMaterial(name),
_conductivity(k)
{
}




fourierMaterial::fourierMaterial(const std::string& name) :
conductorMaterial(name),
_conductivity(0.0)
{
}




bool fourierMaterial::check() const
{
    bool ret = conductorMaterial::check();

    if (_conductivity <= 0.0)
    {
        ret = false;
        std::cout << "Error in fourier material. Non-positive conductivity";
    }

    return ret;
}




conductorMP* fourierMaterial::createMaterialPoint() const
{
    fourierMP *mp = new fourierMP(*this);
    return mp;
}




double fourierMaterial::getProperty(const propertyName p) const
{
    double ret=0.0;

    if      (p == PR_CONDUCTIVITY)        ret = _conductivity;
    else if (p == PR_THERMAL_CAP)         ret = _capacity;
    else
    {
        std:: cout << "\n Error in GetfourierMaterialProperty.";
    }
    return ret;
}




void fourierMaterial::print(std::ostream &of) const
{
    of  << "\n   Isotropic thermally conducting material "
        << "\n   Conductivity      : " << _conductivity
        << "\n   Density           : " << density()
        << "\n   Ref. Heat capacity: " << _capacity;
}




void fourierMaterial::setRandom()
{
    conductorMaterial::setRandom();
    _conductivity = muesli::randomUniform(1.0, 10.0);
    _capacity     = muesli::randomUniform(1.0, 10.0);
}




bool fourierMaterial::test(std::ostream& os)
{
    bool isok = true;
    setRandom();

    conductorMP* p = this->createMaterialPoint();

    isok = p->testImplementation(os);
    delete p;
    return isok;
}




fourierMP::fourierMP(const fourierMaterial& theMaterial)
:
    conductorMP(theMaterial),
    mat(&theMaterial)
{
    gradT_n.setZero();
    gradT_c.setZero();
}




void fourierMP::conductivity(istensor& K) const
{
    K = istensor::identity() * mat->_conductivity;
}




void fourierMP::contractTangent(const ivector& na, const ivector& nb, double& tg) const
{
    tg = mat->_conductivity * na.dot(nb);
}




double fourierMP::density() const
{
    return theConductor->density();
}




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

    state.theTime = time_n;
    state.theVector.push_back(gradT_n);

    return state;
}




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

    state.theTime = time_c;
    state.theVector.push_back(gradT_c);

    return state;
}




double fourierMP::heatCapacity() const
{
    return mat->_capacity;
}




double fourierMP::heatCapacityDerivative() const
{
    return 0.0;
}




void fourierMP::heatflux(ivector &q) const
{
    q = (- mat->_conductivity) * gradT_c;
}




void fourierMP::heatfluxDerivative(ivector &qprime) const
{
    qprime.setZero();
}




void fourierMP::setRandom()
{
    conductorMP::setRandom();
}




double fourierMP::thermalEnergy() const
{
     return 0.5*mat->_conductivity * gradT_c.dot(gradT_c);
}




anisotropicConductorMaterial::anisotropicConductorMaterial(const std::string& name,
                                                           const materialProperties& cl)
:
conductorMaterial(name, cl),
_capacity(0.0)
{
    itensor t;
    muesli::assignValue(cl, "kxx", t(0,0));
    muesli::assignValue(cl, "kxy", t(0,1));
    muesli::assignValue(cl, "kxz", t(0,2));
    muesli::assignValue(cl, "kyx", t(1,0));
    muesli::assignValue(cl, "kyy", t(1,1));
    muesli::assignValue(cl, "kyz", t(1,2));
    muesli::assignValue(cl, "kzx", t(2,0));
    muesli::assignValue(cl, "kzy", t(2,1));
    muesli::assignValue(cl, "kzz", t(2,2));

    _K = istensor::symmetricPartOf(t);
    muesli::assignValue(cl, "capacity", _capacity);
}




anisotropicConductorMaterial::anisotropicConductorMaterial(const std::string& name,
                                                           double kxx, double kxy, double kxz,
                                                           double kyy, double kyz, double kzz,
                                                           double rho, double capacity)
:
conductorMaterial(name),
_capacity(capacity)
{
    _K(0,0) = kxx;
    _K(0,1) = _K(1,0) = kxy;
    _K(0,2) = _K(2,0) = kxz;
    _K(1,1) = kyy;
    _K(1,2) = _K(2,1) = kyz;
    _K(2,2) = kzz;
 }




anisotropicConductorMaterial::anisotropicConductorMaterial(const std::string& name) :
conductorMaterial(name),
_capacity(0.0)
{
    _K.setZero();
}




bool anisotropicConductorMaterial::check() const
{
    bool ret = conductorMaterial::check();

    ivector ev = _K.eigenvalues();
    if ( ev.min() <= 0.0)
    {
        ret = false;
        std::cout << "Error in anisotropicConductor material. Non-positive conductivity";
    }

    return ret;
}




conductorMP* anisotropicConductorMaterial::createMaterialPoint() const
{
    anisotropicConductorMP *mp = new anisotropicConductorMP(*this);
    return mp;
}




double anisotropicConductorMaterial::getProperty(const propertyName p) const
{
    double ret=0.0;
    return ret;
}




void anisotropicConductorMaterial::print(std::ostream &of) const
{
    of  << "\n   Anisotropic thermally conducting material "
    << "\n   Conductivity kxx   : " << _K(0,0)
    << "\n   Conductivity kxy   : " << _K(0,1)
    << "\n   Conductivity kxz   : " << _K(0,2)
    << "\n   Conductivity kyy   : " << _K(1,1)
    << "\n   Conductivity kyz   : " << _K(1,2)
    << "\n   Conductivity kzz   : " << _K(2,2)
    << "\n   Density            : " << density()
    << "\n   Ref. Heat capacity : " << _capacity;
}




void anisotropicConductorMaterial::setRandom()
{
    conductorMaterial::setRandom();
    _capacity = muesli::randomUniform(1.0, 10.0);

    _K.setRandom();
}




bool anisotropicConductorMaterial::test(std::ostream& os)
{
    bool isok = true;
    setRandom();

    conductorMP* p = this->createMaterialPoint();

    isok = p->testImplementation(os);
    delete p;

    return isok;
}




anisotropicConductorMP::anisotropicConductorMP(const anisotropicConductorMaterial& theMaterial)
:
conductorMP(theMaterial),
mat(&theMaterial)
{
    gradT_n.setZero();
    gradT_c.setZero();
}




void anisotropicConductorMP::conductivity(istensor& K) const
{
    K = mat->_K;
}




void anisotropicConductorMP::contractTangent(const ivector& na, const ivector& nb, double& tg) const
{
    tg = na.dot(mat->_K * nb);
}




double anisotropicConductorMP::density() const
{
    return theConductor->density();
}




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

    state.theTime = time_n;
    state.theVector.push_back(gradT_n);

    return state;
}




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

    state.theTime = time_c;
    state.theVector.push_back(gradT_c);

    return state;
}




double anisotropicConductorMP::heatCapacity() const
{
    return mat->_capacity;
}




double anisotropicConductorMP::heatCapacityDerivative() const
{
    return 0.0;
}




void anisotropicConductorMP::heatflux(ivector &q) const
{
    q = - (mat->_K * gradT_c);
}




void anisotropicConductorMP::heatfluxDerivative(ivector &qprime) const
{
    qprime.setZero();
}




void anisotropicConductorMP::setRandom()
{
    conductorMP::setRandom();
}




double anisotropicConductorMP::thermalEnergy() const
{
    return 0.5* gradT_c.dot(mat->_K * gradT_c);
}




nonlinearConductorMaterial::nonlinearConductorMaterial(const std::string& name,
                                                       const materialProperties& cl)
:
conductorMaterial(name, cl)
{
    muesli::assignValue(cl, "k0", k0);
    muesli::assignValue(cl, "k1", k1);
    muesli::assignValue(cl, "ts", tS);
    muesli::assignValue(cl, "tl", tL);
    muesli::assignValue(cl, "capacity", _capacity);
}




nonlinearConductorMaterial::nonlinearConductorMaterial(const std::string& name,
                                                       const double _k0, const double _k1,
                                                       const double rho)
:
conductorMaterial(name)
{
     k0 = _k0;
     k1 = _k1;
}




nonlinearConductorMaterial::nonlinearConductorMaterial(const std::string& name) :
conductorMaterial(name)
{
     k0 = 0.0;
     k1 = 0.0;
}




nonlinearConductorMaterial::nonlinearConductorMaterial(const nonlinearConductorMaterial* nlCM) :
conductorMaterial("none")
{
     k0 = nlCM->k0;
     k1 = nlCM->k1;
}




bool nonlinearConductorMaterial::check() const
{
    bool ret = conductorMaterial::check();

    if (k0 <= 0.0)
    {
        ret = false;
        std::cout << "Error in nonlinear material. Non-positive conductivity at temp=0";
    }
    if (k1 <= 0.0)
    {
        ret = false;
        std::cout << "Error in nonlinear material. Non-positive conductivity";
    }

    return ret;
}




conductorMP* nonlinearConductorMaterial::createMaterialPoint() const
{
    nonlinearConductorMP *mp = new nonlinearConductorMP(*this);
    return mp;
}




double nonlinearConductorMaterial::getProperty(const propertyName p) const
{
    double ret=0.0;

    switch (p)
    {
        case PR_CONDUCTIVITY: ret = k0;      break;
        case PR_THERMAL_CAP:  ret = _capacity; break;
        default: ret = 0.0;
    }

    return ret;
}




void nonlinearConductorMaterial::print(std::ostream &of) const
{
    of  << "\n   Nonlinear thermally conducting material "
        << "\n   Conductivity k0   : " << k0
        << "\n   Conductivity k1   : " << k1
        << "\n   Solid temp tS     : " << tS
        << "\n   Liquid temp tL    : " << tL
        << "\n   Density           : " << density()
        << "\n   Ref. heat capacity: " << _capacity;
}




void nonlinearConductorMaterial::setRandom()
{
    conductorMaterial::setRandom();
    k0 = muesli::randomUniform(1.0, 10.0);
    k1 = muesli::randomUniform(0.1, 0.3);
    tS = muesli::randomUniform(1000.0, 1200.0);
    tL = tS + muesli::randomUniform(100.0, 200.0);
    _capacity = muesli::randomUniform(200.0, 300.0);
}




bool nonlinearConductorMaterial::test(std::ostream& os)
{
    bool isok = true;
    setRandom();

    conductorMP* p = this->createMaterialPoint();

    isok = p->testImplementation(os);
    delete p;
    return isok;
}




nonlinearConductorMP::nonlinearConductorMP(const nonlinearConductorMaterial& theMaterial)
:
    conductorMP(theMaterial),
    mat(&theMaterial),
    isPowder(true),
    isFluid(false),
    isSolid(false),
    isFluids(false),
    isSolidf(false)
{
    gradT_n.setZero();
    gradT_c.setZero();
}




nonlinearConductorMP::nonlinearConductorMP(const nonlinearConductorMP& theMaterial)
:
    conductorMP(theMaterial),
    isPowder(true),
    isFluid(false),
    isSolid(false),
    isFluids(false),
    isSolidf(false)
{
    gradT_n.setZero();
    gradT_c.setZero();
}




void nonlinearConductorMP::conductivity(istensor& K) const
{
    double ktheta = 0.0;
    if (temp_c <= mat->tS)
    {
        ktheta = mat->k0 + mat->k1 * temp_c;
    }
    else if (temp_c <= mat->tL)
    {
        ktheta = 125.43472 - 0.0617 * temp_c;
    }
    else
    {
        ktheta = 4.952 + 0.0136 * temp_c;
    }
    K = istensor::scaledIdentity(ktheta);
}




void nonlinearConductorMP::contractTangent(const ivector& na, const ivector& nb, double& tg) const
{
    istensor K;
    conductivity(K);
    tg = na.dot(K*nb);
}




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

    state.theTime = time_n;
    state.theVector.push_back(gradT_n);

    return state;
}




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

    state.theTime = time_c;
    state.theVector.push_back(gradT_c);

    return state;
}




double nonlinearConductorMP::density() const
{
    double rho = 3191.5;
    if (!isPowder)
        rho = theConductor->density();
    return rho;
}




double nonlinearConductorMP::heatCapacity() const
{
    double c = 0.0;
    if (temp_c <= mat->tS)
    {
        c = 0.18 * temp_c + 391.0;
    }
    else if (temp_c <= mat->tL)
    {
        c = 1.2 * temp_c - 1162.0;
    }
    else
    {
        c = 758.0;
    }
    return c;
}




double nonlinearConductorMP::heatCapacityDerivative() const
{
    double cprime = 0.0;
    if (temp_c <= mat->tS)
    {
        cprime = 0.18;
    }
    else if (temp_c <= mat->tL)
    {
        cprime = 1.2;
    }
    else
    {
        cprime = 0.0;
    }
    return cprime;
}




void nonlinearConductorMP::heatflux(ivector &q) const
{
    istensor K;
    conductivity(K);
    q = -(K*gradT_c);
}




void nonlinearConductorMP::heatfluxDerivative(ivector &qprime) const
{
    double kp;
    if (temp_c <= mat->tS)
    {
        kp = mat->k1;
    }
    else if (temp_c <= mat->tL)
    {
        kp = -0.0617;
    }
    else
    {
        kp = 0.0136;
    }
    qprime = (-kp)*gradT_c;
}




double nonlinearConductorMP::phase() const
{
    double phase = 0.0;

    if      (isPowder) phase = 0.0;
    else if (isFluid)  phase = 1.0;
    else if (isSolid)  phase = 2.0;
    else if (isFluids) phase = 1.0;
    else if (isSolidf) phase = 2.0;

    return phase;
}




void nonlinearConductorMP::setRandom()
{
    conductorMP::setRandom();
}




double nonlinearConductorMP::thermalEnergy() const
{
    istensor K;
    conductivity(K);
    return 0.5 * gradT_c.dot(K*gradT_c);
}




void nonlinearConductorMP::updateCurrentState(double theTime, double temp, const ivector& gradT)
{
    conductorMP::updateCurrentState(theTime, temp, gradT);

    if (isPowder && temp_c > mat->tL)
    {
        isPowder = false;
        isFluid  = true;
    }

    else if (isFluid && temp_c < mat->tS)
    {
        isFluid = false;
        isSolid = true;
    }

    else if (isSolid && temp_c > mat->tL)
    {
        isSolid = false;
        isFluids= true;
    }

    else if (isFluids && temp_c < mat->tS)
    {
        isFluids = false;
        isSolidf = true;
    }

    else if (isSolidf && temp_c > mat->tL)
    {
        isSolidf = false;
        isFluids = true;
    }
}
