
/*
 * Material Model Concrete (Simplified concrete model)
 *
 */

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

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

__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 A         = cmat[6];
    double B         = cmat[7];
    double n         = cmat[8];
    double f_c       = cmat[9];
    double T         = cmat[10];
    double epsf_min  = cmat[11];
    double p_c       = cmat[12];
    double p_l       = cmat[13];
    double epsv_c_l  = cmat[14];
    double D1        = cmat[15];
    double D2        = cmat[16];
    double K         = cmat[17];
    double shear     = cmat[42];
    int erosion_flag = __double2int_rn(cmat[80]);

    double stress[6], dstrain[6];

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

    // Load history device from global memory
    double epsp = data.dp_history[offsetHist + 0];
    double epsv_c = data.dp_history[offsetHist + 1];
    double damage = data.dp_history[offsetHist + 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.dp_strain[offset + 0] + data.dp_strain[offset + 1] + data.dp_strain[offset + 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.dp_stress[offset + i] = (i < 3) ? stress[i] - pressure : stress[i];
    }

    data.dp_history[offsetHist + 0] = epsp + depsp;
    data.dp_history[offsetHist + 1] = epsv_c;
    data.dp_history[offsetHist + 2] = damage;
}

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

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_concrete
