
/*
 * Material Model Rubber
 *
 */

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

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

__constant__ double cmat[200];

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

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

    if (idx >= data.num_tasks) return;

    // Load material parameters from constant memory
    double C1   = cmat[6];
    double C2   = cmat[7];
    double bulk = cmat[43];

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

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

    // Volumetric strain and pressure
    double evol = data.dp_strain[offset + 0] + data.dp_strain[offset + 1] + data.dp_strain[offset + 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.dp_history[offsetHist + i] = eval[i];
    }

    // Save new stress to global memory
    for (int i = 0; i < 6; ++i) {
        data.dp_stress[offset + i] = 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_rubber
