
/*
 * Material Model Orthotropic
 *
 */

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

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

__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 c        = cmat[14];
    double cdec     = cmat[15];
    double Xt       = cmat[16];
    double Xc       = cmat[17];
    double Yt       = cmat[18];
    double Yc       = cmat[19];
    double beta     = cmat[20];
    double S        = cmat[21];
    double residual = cmat[23];
    double c11      = cmat[31];
    double c22      = cmat[32];
    double c12      = cmat[33];
    double c23      = cmat[34];
    double c44      = cmat[35];
    double c55      = cmat[36];
    double c33      = cmat[32];
    double c31      = cmat[33];
    double c66      = cmat[35];

    int erosion_flag = __double2int_rn(cmat[80]);

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

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

    // Load fiber directions from global memory
    for (int i = 0; i < 9; ++i) {
        fiber_dir[i] = data.dp_history[offsetHist + 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 from global memory if viscous stress has been defined
        for (int i = 0; i < 6; ++i) {
            dstrain[i] = (i < 3) ? data.dp_dstrain[offset + i] : data.dp_dstrain[offset + 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.dt / cdec) : 0.0;

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

        // Save bfac for stiffness calculation
        data.dp_bfac[idx] = bfac;
    }

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

    // Load failure parameters from global memory
    for (int i = 0; i < 4; ++i) {
        fail[i] = data.dp_history[offsetHist + 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.dp_internal_fail[idx] = 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.dp_stress[offset + i] = tensor_new[i];
    }

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

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

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

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_orthotropic
