#include "mat_metal.h"
#include "mat_user_functions.h"
#include "kernel_mat_metal.h"

namespace user_metal
{
//-------------------------------------------------------------------------------------------------

MatMetal::MatMetal()
{
    // Enable GPU implementation of this material model
    gpu_enabled = true;

    // Material configuration
    num_hist          =  6;    // Number of state variables per integration point
    pos_epsp          =  0;    // For prop damage, sensor output
    pos_temp          =  4;    // For prop damage, prop thermal
    pos_texture       = -1;    // For texture, not used by this material model
    pos_damage        =  1;    // For damage, position in the state variable array
    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         =  5;    // For prop damage
    pos_rate          = -1;    // For prop damage
    pos_evol          = -1;    // For prop damage
    pos_damage1       =  1;    // For prop damage
    pos_damage2       =  2;    // For prop damage
    pos_damage_active = -1;    // For prop damage

    // Curve definitions
    curve.push_back(6);     // A curve ID is present at position 6 in the material parameter array (cmat)

    // Add the location in the state variable array and the output name for contour plot
    contour.emplace_back(0, "Effective plastic strain");
    contour.emplace_back(1, "Damage");
    contour.emplace_back(4, "Temperature");

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

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

void MatMetal::runInit(UserMatInitCPU data) const
{
    double young = data.p_cmat[1];
    double pr    = data.p_cmat[2];
    double bulk  = young / (3.0 * (1.0 - 2.0 * pr));
    double G     = young / (2.0 * (1.0 + pr));

    // Set material properties
    data.p_cmat[43] = bulk;

    // Stress-strain curve
    int idlc = static_cast<int>(round(data.p_cmat[6]));

    // Yield stress
    double sigy0 = data.p_history[3];

    // Set initial yield stress (if not already defined)
    if (sigy0 == 0.0) {
        // load material properties from database
        sigy0 = mat::load_curve(data.p_curve_data, data.p_curve_val, idlc, 0.0);
        data.p_history[3] = sigy0;
    }

    // Stiffness matrix
    data.p_stiffness[0] = G;
    data.p_stiffness[1] = bulk;
}

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

void MatMetal::runMat(UserMatCPU data) const
{
    // Load material parameters
    double dens  = data.p_cmat[0];
    double epsf  = data.p_cmat[7];
    double shear = data.p_cmat[42];
    double bulk  = data.p_cmat[43];
    double Cp    = data.p_cmat[71];
    double k     = data.p_cmat[73];

    int idlc = static_cast<int>(round(data.p_cmat[6]));

    // Load stress and dstrain data
    double stress[6], dstrain[6];
    for (int i = 0; i < 6; ++i) {
        stress[i] = data.p_stress[i];
        dstrain[i] = data.p_dstrain[i];
    }

    // Load history data
    double epsp        = data.p_history[0];
    double damage1     = data.p_history[1];
    double damage2     = data.p_history[2];
    double sigy0       = data.p_history[3];
    double temperature = data.p_history[4];

    // Remove pressure from old deviatoric stress
    double pressure = -(stress[0] + stress[1] + stress[2]) / 3.0;
    for (int i = 0; i < 3; ++i) {
        stress[i] += pressure;
    }

    // Decomposition of total strain increment into volumetric and deviatoric part
    double deps_vol = (dstrain[0] + dstrain[1] + dstrain[2]) / 3.0;
    for (int i = 0; i < 6; ++i) {
        double dstrain_dev = (i < 3) ? dstrain[i] - deps_vol : dstrain[i] / 2.0;
        stress[i] += 2.0 * shear * dstrain_dev;
    }

    // Pressure based on volumetric strain or EOS
    pressure = -bulk * (data.p_strain[0] + data.p_strain[1] + data.p_strain[2]);

    int yield = 0;
    double deps = 0.0;
    double depsp = 0.0;
    double sig_eff = 0.0;

    // Von Mises yield criterion
    mat::mat_yield_von_mises_1(shear, sigy0, &sig_eff, stress, &deps, &yield);

    if (yield == 1) {
        // Evaluate stress-strain curve
        double sigy1 = mat::load_curve(data.p_curve_data, data.p_curve_val, idlc, epsp + deps);

        // Tangential hardening modulus
        double hardening = (sigy1 - sigy0) / deps;

        // Radial return (returns updated stress tensor and plastic strain increment)
        mat::mat_yield_von_mises_2(shear, sigy0, sigy1, &sig_eff, stress, &deps, &epsp, &depsp);

        // Update yield stress and damage
        sigy0 += hardening * depsp;

        if (epsf != 0.0) {
            damage1 += depsp / epsf;
        }
    }

    // Update temperature if we have plastic flow
    if (depsp > 0.0 && Cp != 0.0) {
        temperature += k * sig_eff * depsp / (Cp * dens);
    }

    // Update history
    data.p_history[0] = epsp;
    data.p_history[1] = damage1;
    data.p_history[2] = damage2;
    data.p_history[3] = sigy0;
    data.p_history[4] = temperature;
    data.p_history[5] = depsp;

    // Add pressure and update stress tensor
    for (int i = 0; i < 6; ++i) {
        data.p_stress[i] = (i < 3) ? stress[i] - pressure : stress[i];
    }
}

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

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

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

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

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

} // namespace mat_metal
