#include "mat_rubber.h"
#include "mat_user_functions.h"
#include "kernel_mat_rubber.h"

namespace user_rubber
{
//-------------------------------------------------------------------------------------------------

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

    // Material configuration
    num_hist          =  3;    // 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        = -1;    // For damage, not used by this material model
    erode_flag        =  0;    // 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, "Eigval 1");
    contour.emplace_back(1, "Eigval 2");
    contour.emplace_back(2, "Eigval 3");

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

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

void MatRubber::runInit(UserMatInitCPU data) const
{
    double dens = data.p_cmat[0];
    double bulk = data.p_cmat[1];
    double C1   = data.p_cmat[6];
    double C2   = data.p_cmat[7];

    double young = 6.0 * (C1 + C2);
    double pr    = 0.5 - (young / bulk) / 4.0;
    double G     = young / (2.0 * (1.0 + pr));

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

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

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

void MatRubber::runMat(UserMatCPU data) const
{
    // Load material parameters
    double C1   = data.p_cmat[6];
    double C2   = data.p_cmat[7];
    double bulk = data.p_cmat[43];

    double fmat[9], c[6], eval[3];

    // Load deformation gradient matrix
    for (int i = 0; i < 9; ++i) {
        fmat[i] = data.p_f_mat[i];
    }

    // Volumetric strain and pressure
    double evol = data.p_strain[0] + data.p_strain[1] + data.p_strain[2];
    double pressure = -bulk * evol;

    mat::cauchy_green_tensor(fmat, c);

    // Re-using fmat array as eigen vector array since fmat array is no longer needed
    double* evec = fmat;

    // Get principal stretches (squared)
    mat::calc_eigen_values(c, eval);
    mat::calc_eigen_vectors(c, eval, evec);

    double alpha = 2.0 * C1 * (eval[0] - eval[2]) - 2.0 * C2 * (1.0 / eval[0] - 1.0 / eval[2]);
    double beta = 2.0 * C1 * (eval[1] - eval[2]) - 2.0 * C2 * (1.0 / eval[1] - 1.0 / eval[2]);

    eval[0] = (2.0 * alpha - 1.0 * beta) / 3.0 - pressure;
    eval[1] = (-1.0 * alpha + 2.0 * beta) / 3.0 - pressure;
    eval[2] = (-1.0 * alpha - 1.0 * beta) / 3.0 - pressure;

    // Re-using c array as local stress array since c array is not longer needed
    double* stress = c;
    for (int i = 0; i < 6; ++i) {
        stress[i] = 0.0;
    }

    for (int i = 0; i < 3; ++i) {
        int vector_offset = i * 3;
        stress[0] += eval[i] * pow(evec[vector_offset + 0], 2);
        stress[1] += eval[i] * pow(evec[vector_offset + 1], 2);
        stress[2] += eval[i] * pow(evec[vector_offset + 2], 2);
        stress[3] += eval[i] * evec[vector_offset + 0] * evec[vector_offset + 1];
        stress[4] += eval[i] * evec[vector_offset + 1] * evec[vector_offset + 2];
        stress[5] += eval[i] * evec[vector_offset + 2] * evec[vector_offset + 0];
    }

    // Save eigen values for contour plot output
    for (int i = 0; i < 3; ++i) {
        data.p_history[i] = eval[i];
    }

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

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

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

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

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

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

} // namespace mat_rubber
