#include "mat_orthotropic.h"
#include "mat_user_functions.h"
#include "kernel_mat_orthotropic.h"

namespace user_orthotropic
{
//-------------------------------------------------------------------------------------------------

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

    // Material configuration
    num_hist          = 20;    // 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       = 10;    // For texture, assign starting position 10 in the state variable array
    pos_damage        = 19;    // For damage, position in the state variable array
    erode_flag        =  1;    // Erosion flag (0=no erosion, 1=element erosion, 2=node splitting)
    pos_erode_flag    = 22;    // Copies the defined position to to position 80 in the material parameter array
    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(6, "Damage TF");
    contour.emplace_back(7, "Damage CF");
    contour.emplace_back(8, "Damage TM");
    contour.emplace_back(9, "Damage CM");

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

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

void MatOrthotropic::runInit(UserMatInitCPU data) const
{
    // Get material parameters from input data
    double dens     = data.p_cmat[0];  // Density of the material
    double E1       = data.p_cmat[6];  // Young's modulus in first direction
    double E2       = data.p_cmat[7];  // Young's modulus in second direction
    double G12      = data.p_cmat[8];  // Shear modulus in plane
    double pr12     = data.p_cmat[9];  // Poisson's ratio in plane
    double pr23     = data.p_cmat[10]; // Poisson's ratio transverse
    double c        = data.p_cmat[14]; // Viscosity
    double cdec     = data.p_cmat[15]; // Viscosirt decay coefficient
    double Xt       = data.p_cmat[16]; // Tensile strength in first direction
    double Xc       = data.p_cmat[17]; // Compressive strength in first direction
    double Yt       = data.p_cmat[18]; // Tensile strength in second direction
    double Yc       = data.p_cmat[19]; // Compressive strength in second direction
    double S        = data.p_cmat[21]; // Parameter controlling shear modulus reduction
    double residual = data.p_cmat[23]; // Residual shear strength

    // Calculate intermediate parameters for material model
    double shear = max(max(E1, E2) / (2.0 * (1.0 + min(pr12, pr23))), G12);
    double bulk = max(E1, E2) / (3.0 * (1.0 - 2.0 * max(pr12, pr23)));
    double xi = c / sqrt(2.0 * E2 * dens);

    // Set default values for material parameters if not provided
    if (Xt == 0.0) data.p_cmat[16] = 1.0e20;
    if (Xc == 0.0) data.p_cmat[17] = 1.0e20;
    if (Yt == 0.0) data.p_cmat[18] = 1.0e20;
    if (Yc == 0.0) data.p_cmat[19] = 1.0e20;
    if (S  == 0.0) data.p_cmat[21] = 1.0e20;

    double G23 = E2 / (2.0 * (1.0 + pr23));

    // Create a 3x3 matrix for the material model
    double A[9], Ainv[9];
    A[0] = 1.0 / E1;
    A[4] = 1.0 / E2;
    A[8] = 1.0 / E2;
    A[1] = A[3] = -pr12 / E1;
    A[5] = A[7] = -pr23 / E2;
    A[6] = A[2] = -pr12 / E1;

    // Invert the matrix
    mat::invert_3x3(A, Ainv);

    // Copy inverted matrix to material parameters
    data.p_cmat[31] = Ainv[0];
    data.p_cmat[32] = Ainv[4];
    data.p_cmat[33] = Ainv[3];
    data.p_cmat[34] = Ainv[7];
    data.p_cmat[35] = 2.0 * G12;
    data.p_cmat[36] = 2.0 * G23;

    // Set stiffness values
    data.p_stiffness[0] = shear;
    data.p_stiffness[1] = bulk;
    data.p_stiffness[2] = xi;
}

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

void MatOrthotropic::runMat(UserMatCPU data) const
{
    // Load material parameters from constant memory
    double c        = data.p_cmat[14];
    double cdec     = data.p_cmat[15];
    double Xt       = data.p_cmat[16];
    double Xc       = data.p_cmat[17];
    double Yt       = data.p_cmat[18];
    double Yc       = data.p_cmat[19];
    double beta     = data.p_cmat[20];
    double S        = data.p_cmat[21];
    double residual = data.p_cmat[23];
    double c11      = data.p_cmat[31];
    double c22      = data.p_cmat[32];
    double c12      = data.p_cmat[33];
    double c23      = data.p_cmat[34];
    double c44      = data.p_cmat[35];
    double c55      = data.p_cmat[36];
    double c33      = data.p_cmat[32];
    double c31      = data.p_cmat[33];
    double c66      = data.p_cmat[35];

    int erosion_flag = static_cast<int>(round(data.p_cmat[80]));

    double stress[6], strain[6], tensor_new[6], fiber_dir[9];

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

    // Load fiber directions
    for (int i = 0; i < 9; ++i) {
        fiber_dir[i] = data.p_history[10 + i];
    }

    // Transform strain from material to fiber directions
    mat::transform_tensor<double, 1>(strain, fiber_dir, tensor_new);

    stress[0] = c11 * tensor_new[0] + c12 * tensor_new[1] + c31 * tensor_new[2];
    stress[1] = c12 * tensor_new[0] + c22 * tensor_new[1] + c23 * tensor_new[2];
    stress[2] = c31 * tensor_new[0] + c23 * tensor_new[1] + c33 * tensor_new[2];
    stress[3] = c44 * tensor_new[3];
    stress[4] = c55 * tensor_new[4];
    stress[5] = c66 * tensor_new[5];

    // Add viscous stresses
    if (c != 0.0) {
        // Re-using strain array as dstrain as strain is not used anymore
        double* dstrain = strain;

        // Load dstrain
        for (int i = 0; i < 6; ++i) {
            dstrain[i] = (i < 3) ? data.p_dstrain[i] : data.p_dstrain[i] / 2.0;
        }

        // Transform dstrain from material to fiber directions
        mat::transform_tensor<double, 1>(dstrain, fiber_dir, tensor_new);

        double bfac = (cdec != 0.0) ? exp(-data.dt1 / cdec) : 0.0;

        for (int i = 0; i < 6; ++i) {
            double srat = bfac * data.p_history[i] + (1.0 - bfac) * tensor_new[i] / data.dt1;
            stress[i] += c * srat;
            data.p_history[i] = srat;
        }

        // Save bfac for stiffness calculation
        data.p_stiffness[3] = max(bfac, data.p_stiffness[3]);
    }

    // Re-using strain array as fail as strain is not used anymore
    double* fail = strain;

    // Load failure parameters
    for (int i = 0; i < 4; ++i) {
        fail[i] = data.p_history[6 + i];
    }

    // Check failure
    double tmax = fabs(stress[3]);
    double dmg = 0.0;

    // Fiber in tension/shear
    if (stress[0] > 0.0) {
        dmg = min(1.0, pow(stress[0] / Xt, 2) + beta * (tmax / S));
        fail[0] = max(fail[0], dmg);
    }
    // Fiber in compression/shear
    else {
        dmg = min(1.0, pow(stress[0] / Xc, 2));
        fail[1] = max(fail[1], dmg);
    }

    // Matrix tension/shear
    double smax = stress[1];
    if (smax > 0.0) {
        dmg = min(1.0, pow(smax / Yt, 2) + beta * (tmax / S));
        fail[2] = max(fail[2], dmg);
    }

    // Matrix compression/shear
    double smin = stress[1];
    if (smin < 0.0) {
        dmg = pow(smax / (2.0 * S), 2) + (pow(Yc / (2.0 * S), 2) - 1.0) * smax / Yc +
              pow(tmax / S, 2);
        dmg = min(1.0, dmg);
        fail[3] = max(fail[3], dmg);
    }

    // Tensile fiber failure => remove all stresses
    if (fail[0] == 1.0 && erosion_flag <= 1.0) {
        for (int i = 0; i < 6; ++i) {
            stress[i] *= residual;
        }

        if (erosion_flag == 1.0) {
            *data.p_internal_fail += 1;
        }
    }

    // Fiber buckling => remove compressive fiber stresses and in-plane shear stress component
    if (fail[1] == 1.0) {
        stress[0] = max(stress[0], residual * stress[0]);
        stress[3] = residual * stress[3];
    }

    // Matrix failure => remove tensile stresses in matrix and all in-plane shear stress
    if (fail[2] == 1.0 || fail[3] == 1.0) {
        stress[1] = min(stress[1], residual * stress[1]);
        stress[3] = residual * stress[3];
    }

    // Transform stress from fiber to material directions
    mat::transform_tensor<double, 0>(stress, fiber_dir, tensor_new);

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

    for (int i = 0; i < 4; ++i) {
        data.p_history[6 + i] = fail[i];
    }

    data.p_history[19] = max(max(fail[0], fail[1]), max(fail[2], fail[3]));
}

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

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

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

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

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

} // namespace user_orthotropic
