%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-10-29 15:04:15
 % @LastEditTime: 2021-04-19 10:11:20
%--------------------------
% gradient of shape function on one ele, four nodes, N1 ~ N4
%
% gradN: [4*8]
% [dN1dx, 0,     dN2dx, 0, ...
%  0,     dN1dx, 0,     dN2dx, ...
%  dN1dy, 0,     dN2dy, 0, ...
%  0,     dN1dy, 0,     dN2dy, ...]
% This sequence is setted for computaion
%
% gradN2: [3*8]
% [dN1dx, 0,     dN2dx, 0,     ...
%  0,     dN1dy, 0,     dN2dy, ...
%  dN1dy, dN1dx, dN2dy, dN2dx]
function [gradN, gradN2, N] = shape_gradient_ele(a, b, sample_num)
    if sample_num==1
        [gradN] = shape_gradient_1(a, b);
        N = [];
        gradN2 = [];
    else
        [gradN, gradN2, N] = shape_gradient_4(a, b);
    end
end

% 4*8
function [gradN] = shape_gradient_1(a, b)    
    GN_x = 0; 
    GN_y = 0;            
    
    dN_x = 1/4*[-(1-GN_x)  (1-GN_x) (1+GN_x) -(1+GN_x)];            
    dN_y = 1/4*[-(1-GN_y) -(1+GN_y) (1+GN_y)  (1-GN_y)];   

    J = [dN_x; dN_y]*[ -a  a  a  -a;  
                       -b  -b  b  b]';            
    G = [inv(J) zeros(size(J)); 
        zeros(size(J)) inv(J)];      

    dN(1,1:2:8) = dN_x; dN(3,1:2:8) = dN_y;            
    dN(2,2:2:8) = dN_x; dN(4,2:2:8) = dN_y;   

    B = G*dN;  
    gradN = B;
end

% 4 * 8 & 3*8 for every cell
function [gradN, gradN2, N] = shape_gradient_4(a, b)    
    GaussNodes = [-1/sqrt(3); 1/sqrt(3)];
    gradN  = cell(4,1);
    gradN2 = cell(4,1);
    L = [1 0 0 0; 0 0 0 1; 0 1 1 0];

    nodes = [-1,  1];
    N = cell(4,1);
    
    for i = 1:2
        for j = 1:2
            gp = 2*(i-1)+j;       
            GN_x = GaussNodes(i); 
            GN_y = GaussNodes(j);  
            % GN_x = nodes(i); GN_y = nodes(j);          
            
            dN_x = 1/4*[-(1-GN_x)  (1-GN_x) (1+GN_x) -(1+GN_x)];            
            dN_y = 1/4*[-(1-GN_y) -(1+GN_y) (1+GN_y)  (1-GN_y)];   

            J = [dN_x; dN_y]*[ -a  a  a  -a;  
                               -b  -b  b  b]';          
            G = [inv(J) zeros(size(J)); 
                zeros(size(J)) inv(J)];  
            
            dN1(1,1:2:8) = dN_x; dN1(2,1:2:8) = dN_y;            
            dN1(3,2:2:8) = dN_x; dN1(4,2:2:8) = dN_y; 
                                    
            % dNdx = 4*8
            B = G*dN1;  
            gradN{gp} = B([1,3,2,4],:);         

            % dNdx = 3*8
            dN1(1,1:2:8) = dN_x; dN1(2,1:2:8) = dN_y;            
            dN1(3,2:2:8) = dN_x; dN1(4,2:2:8) = dN_y;   
            B2 = L*G*dN1;  
            gradN2{gp} = B2; 

            % N = 2*8 
            GN_x = nodes(i); GN_y = nodes(j);
            N1 = (1-GN_x)*(1-GN_y)/4;
            N2 = (1+GN_x)*(1-GN_y)/4;
            N3 = (1+GN_x)*(1+GN_y)/4;
            N4 = (1-GN_x)*(1+GN_y)/4;
            N{gp} = [N1, 0, N2, 0, N3, 0, N4, 0;
                    0, N1, 0, N2, 0, N3, 0, N4];
            % N2{gp} = [N1, N1; 
            %         N2, N2; 
            %         N3, N3; 
            %         N4, N4;]; % 4*2
        end
    end
end