
/*
 * Material Model Metal: \m/
 *
 */

#include "gpu_error.h"
#include "kernel_mat_metal.h"
#include "mat_user_functions.h"
#include "mat_user_defines.h"

namespace user_metal
{

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

__constant__ double cmat[200];

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

__global__ void kernel_material(UserMatDevice data)
{
    const int idx = threadIdx.x + blockDim.x * blockIdx.x;
    const int elem_id = idx / data.num_ip;
    const int offset = idx * 6;
    const int offsetHist = idx * data.num_history;

    if (idx >= data.num_tasks) return;
    if (data.dp_eroded[elem_id] > 0) return;

    // Load material parameters from constant memory
    double dens  = cmat[0];
    double epsf  = cmat[7];
    double shear = cmat[42];
    double bulk  = cmat[43];
    double Cp    = cmat[71];
    double k     = cmat[73];

    // Subract with 1, because the defined curve indices starts at index 1
    int idlc = __double2int_rn(cmat[6]);

    // Load stress and dstrain data from global memory
    double stress[6], dstrain[6];
    for (int i = 0; i < 6; ++i) {
        stress[i] = data.dp_stress[offset + i];
        dstrain[i] = data.dp_dstrain[offset + i];
    }

    // Load history data from global memory
    double epsp = data.dp_history[offsetHist + 0];
    double damage1 = data.dp_history[offsetHist + 1];
    double damage2 = data.dp_history[offsetHist + 2];
    double sigy0 = data.dp_history[offsetHist + 3];
    double temperature = data.dp_history[offsetHist + 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 = (data.eos == 0) ? -bulk * (data.dp_strain[offset + 0] + data.dp_strain[offset + 1] + data.dp_strain[offset + 2])
                               : data.dp_eos_pressure[idx];

    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.dp_curve_data, data.dp_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.dp_history[offsetHist + 0] = epsp;
    data.dp_history[offsetHist + 1] = damage1;
    data.dp_history[offsetHist + 2] = damage2;
    data.dp_history[offsetHist + 3] = sigy0;
    data.dp_history[offsetHist + 4] = temperature;
    data.dp_history[offsetHist + 5] = depsp;

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

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

void kernelMaterial(UserMatHost host, UserMatDevice device, cudaStream_t stream)
{
    cudaMemcpyToSymbolAsync(cmat, host.p_cmat, sizeof(double) * 200, 0, cudaMemcpyHostToDevice, stream);

    const unsigned int block_size = 128;
    const unsigned int num_blocks = calcNumBlocks(device.num_tasks, block_size);

    kernel_material<<<num_blocks, block_size, 0, stream>>>(device);
    kernel_error_check(stream, __FILE__, __LINE__);
}

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

} // namespace user_metal
