/****************************************************************************
*
*                                 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 "brownmiller.h"
#include <iostream>
#include <cmath>

using namespace muesli;

brownMillerCriterium::brownMillerCriterium()
:
    E(0.0),
    nu(0.0),
    sigmaf(0.0),
    epsf(0.0),
    alpha(0.0),
b(0.0),
c(0.0)
{

}




brownMillerCriterium::brownMillerCriterium(const materialProperties& cl)
:
    E(0.0),
    nu(0.0),
    sigmaf(0.0),
    epsf(0.0),
    alpha(0.0),
    b(0.0),
    c(0.0)
{
    muesli::assignValue(cl, "young",  E);
    muesli::assignValue(cl, "poisson", nu);
    muesli::assignValue(cl, "sigmaf", sigmaf);
    muesli::assignValue(cl, "epsf",   epsf);
    muesli::assignValue(cl, "alpha",  alpha);
    muesli::assignValue(cl, "b",  b);
    muesli::assignValue(cl, "c",  c);
}




brownMillerCriterium::~brownMillerCriterium()
{

}




bool brownMillerCriterium::check() const
{
    bool ok = true;

    if (E <= 0.0) ok = false;
    if (sigmaf <= 0.0) ok = false;
    if (epsf <= 0.0) ok = false;
    if (alpha < 0.0) ok = false;
    if (alpha > 1.0) ok = false;

    return ok;
}




ssFailurePoint* brownMillerCriterium::createFailurePoint() const
{
    return new brownMillerPoint(*this);
}




double brownMillerCriterium::getProperty(const propertyName p) const
{
    return 0.0;
}




void brownMillerCriterium::print(std::ostream &of) const
{
    of  << "\n\n Kandil-Brown-Miller criterium for critical plane fatigue analysis"
        << "\n   Young modulus              : " << E
        << "\n   Poisson's ratio            : " << nu
        << "\n   Fatigue stress limit sigmaf: " << sigmaf
        << "\n   Fatigue stress exponent b  : " << b
        << "\n   Fatigue strain limit epsf  : " << epsf
        << "\n   Fatigue strain exponent c  : " << c
        << "\n   KBM averaging parameter a  : " << alpha;
}




void brownMillerCriterium::setRandom()
{
    E      = muesli::randomUniform(10.0, 10000.0);
    sigmaf = muesli::randomUniform(1.0, 100.0 );
    epsf   = muesli::randomUniform(1.0e-2, 1e-1);
    alpha  = muesli::randomUniform(0.0, 1.0);
}




brownMillerPoint::brownMillerPoint(const brownMillerCriterium &m)
:
theCriterium(m),
maxGamma_n(-__DBL_MAX__),
minGamma_n(__DBL_MAX__),
maxEpsm_n(-__DBL_MAX__),
minEpsm_n(__DBL_MAX__),
maxGamma_c(-__DBL_MAX__),
minGamma_c(__DBL_MAX__),
maxEpsm_c(-__DBL_MAX__),
minEpsm_c(__DBL_MAX__)
{

}




materialState brownMillerPoint::getConvergedState() const
{
    materialState st;
    st.theDouble.push_back(maxGamma_n - minGamma_n);
    st.theDouble.push_back(maxEpsm_n - minEpsm_n);

    return st;
}




materialState brownMillerPoint::getCurrentState() const
{
    materialState st;
    st.theDouble.push_back(maxGamma_c - minGamma_c);
    st.theDouble.push_back(maxEpsm_c - minEpsm_c);

    return st;
}




void brownMillerPoint::commitCurrentState()
{
    maxGamma_n = std::max<double>(maxGamma_c, maxGamma_n);
    minGamma_n = std::min<double>(minGamma_c, minGamma_c);
    maxEpsm_n  = std::max<double>(maxEpsm_c, maxEpsm_n);
    minEpsm_n  = std::min<double>(minEpsm_c, minEpsm_n);
}




void brownMillerPoint::resetCurrentState()
{
    maxGamma_c = maxGamma_n;
    minGamma_c = minGamma_n;
    maxEpsm_c  = maxEpsm_n;
    minEpsm_c  = minEpsm_n;
}




void brownMillerPoint::updateCurrentState(double t, const istensor& strain,
                                            const istensor& stress)
{
    ivector eig = strain.eigenvalues();

    double epsI = eig(0);   // the smallest
    double epsIII = eig(2); // the largest

    double gamma_max = epsIII - epsI;
    double epsm = 0.5*(epsI + epsIII);

    maxGamma_c = std::max<double>(maxGamma_c, gamma_max);
    minGamma_c = std::min<double>(minGamma_c, gamma_max);

    maxEpsm_c  = std::max<double>(maxEpsm_c, epsm);
    minEpsm_c = std::min<double>(minEpsm_c, epsm);
}




bool brownMillerPoint::isIntact() const
{
    return true;
}




double brownMillerPoint::remainingLife() const
{
    // limits for lime (alog(n))
    double inf = 1.0;
    double sup = 9.0;

    const double eqRange = (maxGamma_c-minGamma_c) + theCriterium.alpha * (maxEpsm_c - minEpsm_c);
    
    const double nu_e  = theCriterium.nu;
    const double nu_p  = 0.5;
    const double alpha = theCriterium.alpha;
    const double sigmaf= theCriterium.sigmaf;
    const double epsf  = theCriterium.epsf;
    const double E     = theCriterium.E;
    const double b     = theCriterium.b;
    const double c     = theCriterium.c;

    const double A     = 1.0 + nu_e + alpha*(1.0 - nu_e);
    const double B     = 1.0 + nu_p + alpha*(1.0 - nu_p);

    auto kbm = [&] (double m)
    {
        double n = pow(10.0, m);
        return A*sigmaf/E* pow(2.0*n, b) + B*epsf*pow(2.0*n,c) - eqRange;
    };

    const double tol = 1e-4;
    double m = 0.5*(inf+sup);

    double y1 = kbm(inf);
    double y2 = kbm(sup);

    // y2+eps_eq is the eps range that the point resists at 10^sup cycles
    // if positive, it means that the eps_eq is smaller than the fatigue limit,
    // and thus it lifes forever at this level
    if (y2 > 0.0) return sup;

    // y1+eps_eq is the eps range that the point resists at 10^inf cicles
    // if y1 is negative, it means that the eps_eq is larger than the admissible strain
    // at the inf number of cycles, thus it breaks
    if (y1 < 0.0) return 0.0;

    if (y1*y2 > 0.0)
    {
        std::cout << "\n Error in KBM. Solution is not within interval";
    }

    // bisection to determine life n
    while (fabs(sup-inf) > tol)
    {
        m = 0.5*(inf+sup);
        double y3 = kbm(m);

        if (y1*y3 > 0.0)
        {
            inf = m;
            y1  = y3;
        }
        else
        {
            sup = m;
        }
    }

    // returns log10 of number of cycles
    return m;
}
