         float J23 = __powf(J, -(float)2/3);   // J23 = J^(-2/3)
         float x1 = J23*Mu;
         // Bracketed term is I4 = A:C
-        float x2 = J23*(A_gpu[0]*C11+A_gpu[1]*C22+A_gpu[2]*C33+2*A_gpu[3]*C12+2*A_gpu[4]*C23+2*A_gpu[5]*C13) - 1;
-        float x3 = J23*Eta_gpu*x2;
-        float x4 = __fdividef(-(Eta_gpu*x2*(x2+1)+ x1*(C11+C22+C33)), 3.0f);
+        float x2 = J23*(A_hex_gpu[0]*C11+A_hex_gpu[1]*C22+A_hex_gpu[2]*C33+2*A_hex_gpu[3]*C12+2*A_hex_gpu[4]*C23+2*A_hex_gpu[5]*C13) - 1;
+        float x3 = J23*Eta_hex_gpu*x2;
+        float x4 = __fdividef(-(Eta_hex_gpu*x2*(x2+1)+ x1*(C11+C22+C33)), 3.0f);
         float K = Lambda + __fdividef(2*Mu, 3.0f);
         float x5 = K*J*(J-1);
 
         /// Elastic component of the response (isochoric part + volumetric part)
         float SiE11, SiE12, SiE13, SiE22, SiE23, SiE33;
-        SiE11 = x3*A_gpu[0] + x4*Ci11 + x1;
-        SiE22 = x3*A_gpu[1] + x4*Ci22 + x1;
-        SiE33 = x3*A_gpu[2] + x4*Ci33 + x1;
-        SiE12 = x3*A_gpu[3] + x4*Ci12;
-        SiE23 = x3*A_gpu[4] + x4*Ci23;
-        SiE13 = x3*A_gpu[5] + x4*Ci13;
+        SiE11 = x3*A_hex_gpu[0] + x4*Ci11 + x1;
+        SiE22 = x3*A_hex_gpu[1] + x4*Ci22 + x1;
+        SiE33 = x3*A_hex_gpu[2] + x4*Ci33 + x1;
+        SiE12 = x3*A_hex_gpu[3] + x4*Ci12;
+        SiE23 = x3*A_hex_gpu[4] + x4*Ci23;
+        SiE13 = x3*A_hex_gpu[5] + x4*Ci13;
 
         float SvE11, SvE12, SvE13, SvE22, SvE23, SvE33;
         SvE11 = x5*Ci11;

 */
void ClearGPU_Visco(void)
{
    mycudaFree(Di1_gpu);
    mycudaFree(Di2_gpu);
    mycudaFree(Dv1_gpu);
    mycudaFree(Dv2_gpu);

    myprintf("Memory on GPU for viscoelasticity cleared\n");
}

#ifdef ISOCHORIC
        // Isochoric part
        temp = tex1Dfetch(Di1_ref, index);
        temp.x *= Ai_hex_gpu[1]; temp.x += Ai_hex_gpu[0]*SiE11;
        SPK[0] -= temp.x;
        temp.y *= Ai_hex_gpu[1]; temp.y += Ai_hex_gpu[0]*SiE22;
        SPK[1] -= temp.y;
        temp.z *= Ai_hex_gpu[1]; temp.z += Ai_hex_gpu[0]*SiE33;
        SPK[2] -= temp.z;
        temp.w *= Ai_hex_gpu[+1]; temp.w += Ai_hex_gpu[0]*SiE12;
        SPK[3] -= temp.w;
        Di1[index] = make_float4(temp.x, temp.y, temp.z, temp.w);

        temp = tex1Dfetch(Di2_ref, index);
        temp.x *= Ai_hex_gpu[1]; temp.x += Ai_hex_gpu[0]*SiE23;
        SPK[4] -= temp.x;
        temp.y *= Ai_hex_gpu[1]; temp.y += Ai_hex_gpu[0]*SiE13;
        SPK[5] -= temp.y;
        Di2[index] = make_float4(temp.x, temp.y, 0, 0);

#else
        // Volumetric part
        temp = tex1Dfetch(Dv1_ref, index);
        temp.x *= Av_hex_gpu[1]; temp.x += Av_hex_gpu[0]*SvE11;
        SPK[0] -= temp.x;
        temp.y *= Av_hex_gpu[1]; temp.y += Av_hex_gpu[0]*SvE22;
        SPK[1] -= temp.y;
        temp.z *= Av_hex_gpu[1]; temp.z += Av_hex_gpu[0]*SvE33;
        SPK[2] -= temp.z;
        temp.w *= Av_hex_gpu[1]; temp.w += Av_hex_gpu[0]*SvE12;
        SPK[3] -= temp.w;
        Dv1[tid] = make_float4(temp.x, temp.y, temp.z, temp.w);

        temp = tex1Dfetch(Dv2_ref, index);
        temp.x *= Av_hex_gpu[1]; temp.x += Av_hex_gpu[0]*SvE23;
        SPK[4] -= temp.x;
        temp.y *= Av_hex_gpu[1]; temp.y += Av_hex_gpu[0]*SvE13;
        SPK[5] -= temp.y;
        Dv2[tid] = make_float4(temp.x, temp.y, 0, 0);
#endif
