#include "mat_concrete.h"
#include "mat_user_functions.h"
#include "kernel_mat_concrete.h"

namespace user_concrete
{
//-------------------------------------------------------------------------------------------------

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

    // Material configuration
    num_hist          =  3;    // Number of state variables per integration point
    pos_epsp          =  0;    // 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        =  2;    // 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         = -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, "Effective plastic strain");
    contour.emplace_back(1, "Inelastic compaction");
    contour.emplace_back(2, "Damage");

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

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

void MatConcrete::runInit(UserMatInitCPU data) const
{
    double G     = data.p_cmat[1];
    double bulk  = data.p_cmat[17];

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

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

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

void MatConcrete::runMat(UserMatCPU data) const
{
    // Load material parameters
    double A         = data.p_cmat[6];
    double B         = data.p_cmat[7];
    double n         = data.p_cmat[8];
    double f_c       = data.p_cmat[9];
    double T         = data.p_cmat[10];
    double epsf_min  = data.p_cmat[11];
    double p_c       = data.p_cmat[12];
    double p_l       = data.p_cmat[13];
    double epsv_c_l  = data.p_cmat[14];
    double D1        = data.p_cmat[15];
    double D2        = data.p_cmat[16];
    double K         = data.p_cmat[17];
    double shear     = data.p_cmat[42];
    int erosion_flag = static_cast<int>(round(data.p_cmat[80]));

    double stress[6], dstrain[6];

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

    // Load history device
    double epsp   = data.p_history[0];
    double epsv_c = data.p_history[1];
    double damage = data.p_history[2];

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

    // Deviatoric strain increment
    double deps_vol = (dstrain[0] + dstrain[1] + dstrain[2]);
    for (int i = 0; i < 6; ++i) {
        double dstrain_dev = (i < 3) ? dstrain[i] - deps_vol / 3.0 : dstrain[i] / 2.0;
        stress[i] += 2.0 * shear * dstrain_dev;
    }

    // Volumetric update
    double epsv = -(data.p_strain[ 0] + data.p_strain[ 1] + data.p_strain[ 2]);

    // Initialize incremental inelastic compaction
    double depsv_c = 0.0;

    // Inelastic compaction
    double xi = max(0.0, epsv_c / epsv_c_l);
    if (xi < 1.0) {
        // Update pressure
        pressure = K * (epsv - epsv_c);

        // Compaction pressure
        double p_flow = (1.0 - xi) * p_c + xi * p_l;

        if (pressure > p_flow) {
            depsv_c = (pressure - p_flow) / K;

            // Check cap
            if ((epsv_c + depsv_c) > epsv_c_l) {
                epsv_c = epsv_c_l;
                depsv_c = epsv_c_l - epsv_c;
                pressure = p_l - K * (deps_vol + depsv_c);
            } else {
                epsv_c += depsv_c;
                pressure = p_flow;
            }
        }
    }
    // Fully compacted
    else {
        pressure = K * (epsv - epsv_c_l);
    }

    // Tensile cut-off
    double p_cut = -(1.0 - damage) * T;
    if (pressure < p_cut) {
        pressure = p_cut;
    }

    // Effective trial stress
    double sig_eff = mat::mat_effective_stress(stress);

    // Normalized pressure and tensile cut-off
    double p_s = pressure / f_c;
    double T_s = T / f_c;

    // Yield strength
    double sigy = (p_s > 0.0) ? f_c * (A * (1.0 - damage) + B * pow(p_s, n))
                              : f_c * A * (1.0 + pressure / T) * (1.0 - damage);

    // Deviatoric plastic flow
    double depsp = 0.0;

    if (sig_eff > sigy) {
        depsp = (sig_eff - sigy) / (3.0 * shear);
        for (int i = 0; i < 6; ++i) {
            stress[i] *= sigy / sig_eff;
        }
    }

    // Failure strain
    double epsf = max(epsf_min, D1 * pow(p_s + T_s, D2));

    // Update damage
    if (depsv_c > 0.0 || depsp > 0.0) {
        if (damage < 1.0) {
            damage = min(1.0, damage + (depsv_c + depsp) / epsf);
        }
    }

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

    // Update history
    data.p_history[0] = epsp + depsp;
    data.p_history[1] = epsv_c;
    data.p_history[2] = damage;
}

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

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

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

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

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

} // namespace mat_concrete
