#include "mat_glass.h"
#include "mat_user_functions.h"

namespace user_glass
{
//-------------------------------------------------------------------------------------------------

MatGlass::MatGlass()
{
    // Enable GPU implementation of this material model
    gpu_enabled = false;

    // Material configuration
    num_hist          = 15;    // Number of state variables per integration point
    pos_epsp          = -1;    // For prop damage, sensor output
    pos_temp          = -1;    // For prop damage, prop thermal
    pos_texture       = -1;    // For texture, not used by this material model
    pos_damage        =  0;    // For damage, not used by this material model
    erode_flag        =  2;    // Erosion flag (0=no erosion, 1=element erosion, 2=node splitting)
    pos_erode_flag    = -1;    // Not used, since this material model does not use element erosion
    pos_depsp         = -1;    // For prop damage
    pos_rate          = -1;    // For prop damage
    pos_evol          = -1;    // For prop damage
    pos_damage1       = -1;    // For prop damage
    pos_damage2       = -1;    // For prop damage
    pos_damage_active = -1;    // For prop damage

    // Curve definitions
    // curve.push_back(X);

    // Add the location in the state variable array and the output name for contour plot
    contour.emplace_back(0, "Damage");

    // Stiffness reduction (for GPU only)
    stiff_shear = NO_REDUCTION;
    stiff_bulk  = NO_REDUCTION;
    stiff_xi    = NO_REDUCTION;
    stiff_bfac  = MIN_REDUCTION;
}

//-------------------------------------------------------------------------------------------------

void MatGlass::runInit(UserMatInitCPU data) const
{
    double dens  = data.p_cmat[0];
    double young = data.p_cmat[1];
    double c     = data.p_cmat[6];
    int flaw_dim = static_cast<int>(data.p_cmat[14]);
    double shear = data.p_cmat[42];
    double bulk  = data.p_cmat[43];

    // Get defined top and bottom surface normals
    double surface_t[3] = { data.p_cmat[15], data.p_cmat[16], data.p_cmat[17] };
    double surface_b[3] = { data.p_cmat[18], data.p_cmat[19], data.p_cmat[20] };

    double flaw     = data.p_history[7];
    double ac_ratio = data.p_history[11];

    double y_max = 0.0;

    if (flaw != 0.0) {
        for (int i = 0; i < 90; ++i) {
            double phi = M_PI * static_cast<float>(i) / 180.0;
            double lambda = (1.13 - 0.09 * ac_ratio) * (1.0 + 0.1 * pow(1.0 - sin(phi), 2));
            double Q = 1.0 + 1.464 * pow(ac_ratio, 1.65);
            double f = pow((pow(sin(phi), 2) + pow(ac_ratio, 2) * pow(cos(phi), 2)), 1.0 / 4.0);
            double y = lambda * f / sqrt(Q);
            y_max = max(y, y_max);
        }
    }

    data.p_history[11] = y_max;


    // Set stiffness values
    data.p_stiffness[0] = shear;
    data.p_stiffness[1] = bulk;
    data.p_stiffness[2] = c / sqrt(2.0 * young * dens);
}

//-------------------------------------------------------------------------------------------------

void MatGlass::runMat(UserMatCPU data) const
{
    double dens     = data.p_cmat[0];
    double c        = data.p_cmat[6];
    double cdec     = data.p_cmat[7];
    double stress_c = data.p_cmat[8];
    double k_c      = data.p_cmat[9];
    double t_c      = data.p_cmat[10];
    double a_c      = data.p_cmat[11];
    double shear    = data.p_cmat[42];
    double bulk     = data.p_cmat[43];
    double bfac     = (cdec != 0) ? exp(-data.dt1 / cdec) : 0.0;

    double strain[6], stress[6], eval[3];

    // Load data
    for (int i = 0; i < 6; ++i) {
        strain[i] = data.p_strain[i];
    }

    // Load history
    double damage = data.p_history[0];
    double flaw   = data.p_history[7];

    // Deviatoric stress
    double evol = strain[0] + strain[1] + strain[2];
    double pressure = -bulk * evol;

    for (int i = 0; i < 6; ++i) {
        if (i < 3) {
            double eps_dev = strain[i] - evol / 3.0;
            stress[i] = 2.0 * shear * eps_dev - pressure;
        } else {
            double eps_dev = strain[i] / 2.0;
            stress[i] = 2.0 * shear * eps_dev;
        }
    }

    // Damage growth
    mat::calc_eigen_values(stress, eval);
    double eval_max = max(eval[0], eval[1], eval[2]);

    if (eval_max > stress_c) {
        damage += (t_c != 0.0) ? (data.dt1 / t_c) * pow(eval_max / stress_c, a_c) : 1.0;
    }

    // Damage growth at flaw
    if (flaw > 0.0) {
        double dir1 = data.p_history[8];
        double dir2 = data.p_history[9];
        double dir3 = data.p_history[10];
        double form = data.p_history[11];

        double t1 = stress[0] * dir1 + stress[3] * dir2 + stress[5] * dir3;
        double t2 = stress[3] * dir1 + stress[1] * dir2 + stress[4] * dir3;
        double t3 = stress[5] * dir1 + stress[4] * dir2 + stress[2] * dir3;

        double stress_n = t1 * dir1 + t2 * dir2 + t3 * dir3;
        double stress_intensity = max(form * stress_n * sqrt(flaw * M_PI), 0.0);

        if (stress_intensity > k_c) {
            damage += (t_c != 0.0) ? (data.dt1 / t_c) * pow(stress_intensity / k_c, a_c) : 1.0;

            data.p_history[12] = stress_intensity;
        }
    }

    if (c > 0.0) {
        double rate[6];
        for (int i = 0; i < 6; ++i) {
            rate[i] = data.p_dstrain[i] / data.dt1;
        }

        double rv = (rate[0] + rate[1] + rate[2]) / 3.0;
        for (int i = 0; i < 3; ++i) {
            rate[i] -= rv;
        }

        for (int i = 0; i < 6; ++i) {
            double stress_v = data.p_history[1 + i];
            stress_v = bfac * stress_v + 2.0 * c * (1.0 - bfac) * rate[i];

            stress[i] += stress_v;
            data.p_history[1 + i] = stress_v;
        }
    }

    for (int i = 0; i < 6; ++i) {
        data.p_stress[i] = stress[i];
    }

    data.p_history[0] = damage;
    data.p_stiffness[3] = bfac;

    double surface_t[3] = { data.p_cmat[15], data.p_cmat[16], data.p_cmat[17] };
    double surface_b[3] = { data.p_cmat[18], data.p_cmat[19], data.p_cmat[20] };
}

//-------------------------------------------------------------------------------------------------

void MatGlass::runInitGPU(UserMatHost host, UserMatDevice device, cudaStream_t stream) const
{
    // Currently not used, reserved for future implementation
}

//-------------------------------------------------------------------------------------------------

void MatGlass::runMatGPU(UserMatHost host, UserMatDevice device, cudaStream_t stream) const
{
    //kernelMaterial(host, device, stream);
}

//-------------------------------------------------------------------------------------------------

} // namespace mat_glass
