%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-02-19 21:19:39
 % @LastEditTime: 2021-02-19 21:52:22

 % gradient of shape function on one ele, 8 nodes, N1 ~ N8
 % dN: [9*24]
 % [dN1dx, 0,     0,     ...
 %  0,     dN1dx, 0,     ...
 %  0,     0,     dN1dx, ...
 %  dN1dy, 0,     0,     ...
 %  0,     dN1dy, 0,     ...
 %  0,     0,     dN1dy, ...
 %  dN1dz, 0,     0,     ...
 %  0,     dN1dz, 0,     ...
 %  0,     0,     dN1dz, ...
 % This sequence is setted for computaion
%--------------------------
function B = shape_grad(a, b, c)
    GN_x=[-1/sqrt(3),1/sqrt(3)]; GN_y=GN_x; GN_z=GN_x; 
    B = cell(2,2,2);

    % L = zeros(6,9);
    % L(1,1) = 1; L(2,5) = 1; L(3,9) = 1;
    % L(4,2) = 1; L(4,4) = 1; L(5,6) = 1;
    % L(5,8) = 1; L(6,3) = 1; L(6,7) = 1;
    
    for i=1:length(GN_x)
        for j=1:length(GN_y)
            for k=1:length(GN_z)
                x = GN_x(i);
                y = GN_y(j);
                z = GN_z(k);
                dNx = 1/8*[-(1-y)*(1-z)  (1-y)*(1-z)  (1+y)*(1-z) -(1+y)*(1-z) -(1-y)*(1+z)  (1-y)*(1+z)  (1+y)*(1+z) -(1+y)*(1+z)];
                dNy = 1/8*[-(1-x)*(1-z) -(1+x)*(1-z)  (1+x)*(1-z)  (1-x)*(1-z) -(1-x)*(1+z) -(1+x)*(1+z)  (1+x)*(1+z)  (1-x)*(1+z)];
                dNz = 1/8*[-(1-x)*(1-y) -(1+x)*(1-y) -(1+x)*(1+y) -(1-x)*(1+y)  (1-x)*(1-y)  (1+x)*(1-y)  (1+x)*(1+y)  (1-x)*(1+y)];
                J = [dNx;dNy;dNz]*[ -a  a  a  -a  -a  a  a  -a ;  
                                    -b  -b  b  b  -b  -b  b  b; 
                                    -c -c -c -c  c  c  c  c]';
                G = [inv(J) zeros(3) zeros(3);
                    zeros(3) inv(J) zeros(3);
                    zeros(3) zeros(3) inv(J)];

                dN(1,1:3:24) = dNx; dN(2,1:3:24) = dNy; dN(3,1:3:24) = dNz;
                dN(4,2:3:24) = dNx; dN(5,2:3:24) = dNy; dN(6,2:3:24) = dNz;
                dN(7,3:3:24) = dNx; dN(8,3:3:24) = dNy; dN(9,3:3:24) = dNz;
                Be = G*dN;
                B{i,j,k} = Be([1,4,7,2,5,8,3,6,9], :);
            end
        end
    end
end