clear; clearvars -global; clc; close all;
restoredefaultpath;
rng('default');
format short

Ev = 20;
Eh = 40;
nuhh = 0.2;
nuvh = 0.25;
Gvh = 15;

% Ev = 20;
% Eh = 20;
% nuhh = 0.2;
% nuvh = 0.2;
% Gvh = Eh/2/(1+nuhh);
% IMPORTANT RELATION!!! nuvh/Ev = nuhv/Eh

S11 = [1/Eh, -nuhh/Eh, -nuvh/Ev; ...
    -nuhh/Eh, 1/Eh, -nuvh/Ev; ...
    -nuvh/Ev, -nuvh/Ev, 1/Ev];
S = [S11,zeros(3,3);zeros(3,3),diag([2*(1+nuhh)/Eh, 1/Gvh, 1/Gvh])]; % Compliance matrix 6*6 (VTI)
C = S\eye(6); % Stiffness matrix 6*6 (VTI)

lambda = C(1,2); mu_T = C(4,4); mu_L = C(5,5);
alpha = C(1,3)-C(1,2); beta = C(3,3) - lambda - 4*mu_L + 2*mu_T - 2*alpha; % Same unit as Ev/Eh/Gvh for these 5 constants

I2 = eye(3); % Second order identity tensor
I4 = zeros(3,3,3,3); % Fourth order identity tensor
Ce = zeros(3,3,3,3); % Stiffness 4th-order tensor (general case, not necessarily VTI)

n = [sqrt(1/6); sqrt(2/6); sqrt(3/6)];

m = n*n';
for i = 1:3
    for j = 1:3
        for k = 1:3
            for l = 1:3
                I4(i, j, k, l) = (i==k)*(j==l);
                Ce(i, j, k, l) = lambda*(i==j)*(k==l) + 2*mu_T*I4(i, j, k, l) ...
                    + alpha*((i==j)*m(k, l) + m(i, j)*(k==l)) + beta*(m(i, j)*m(k, l)) ...
                    + 2*(mu_L - mu_T)*(m(i, k)*(j==l) + (i==k)*m(j, l));
            end
        end
    end
end
clearvars i j k l

% Define a rotation
Q = rotz(180); % rotx(deg), roty(deg)
e = randn(3,3); e = (e + e')/2;
RHS = Q*double_dot(Ce, e)*Q';
LHS = double_dot(Ce, Q*e*Q');

disp(stiffness_to_mat6by6(Ce));

Ce2 = zeros(3,3,3,3);
for i = 1:3
    for j = 1:3
        for k = 1:3
            for l = 1:3
                for a = 1:3
                    for b = 1:3
                        for c = 1:3
                            for d = 1:3
                                Ce2(i,j,k,l) = Ce2(i,j,k,l) + Q(i,a)*Q(j,b)*Q(k,c)*Q(l,d)*Ce(a,b,c,d);
                            end
                        end
                    end
                end
            end
        end
    end
end
clearvars i j k l a b c d

disp(stiffness_to_mat6by6(Ce2) - stiffness_to_mat6by6(Ce));

E = zeros(9,6); E(1,1) = 1; E(5, 2) = 1; E(9, 3) = 1; E(2, 4) = 0.5; E(4, 4) = 0.5; E(3, 5) = 0.5; E(7, 5) = 0.5; E(6, 6) = 0.5; E(8, 6) = 0.5;
index = [1, 5, 9, 4, 7, 8];

DELAS = reshape(Ce, [9,9])*E; % The use of reshape here is correct! (Specific transformation rule (instead of the Voigt notation, use full matrix/vector version)!)
DELAS = DELAS(index, :);  % 6*6 matrix
disp(DELAS - stiffness_to_mat6by6(Ce));

% To summarize, to reduce the number of independent elastic constants, we
% put constraints on the fourth-order tensor Ce, such that Ce2 = Ce under
% certain Q matrices.

LHS2 = double_dot(Ce2, Q*e*Q'); % This is ALWAYS equal to RHS! (No symmetry constraint on Ce)
% In other words, this is just the tensorial transformation law!