function [Lmat]=Lmatrix

% gauss point location
pos=sqrt(3)/3;
% pos=1;
%Gauss point matrix location
P_g=[-pos,-pos,-pos;
    pos,-pos,-pos;
    pos,pos,-pos;
    -pos,pos,-pos;
    -pos,-pos,pos;
    pos,-pos,pos;
    pos,pos,pos;
    -pos,pos,pos];

% gauss point location
pos=1;
% pos=sqrt(3)/3;
%Gauss point matrix location
coord=[-pos,-pos,-pos;
    pos,-pos,-pos;
    pos,pos,-pos;
    -pos,pos,-pos;
    -pos,-pos,pos;
    pos,-pos,pos;
    pos,pos,pos;
    -pos,pos,pos];

% nodal coordinates, unchanged??
xi_node=coord(:,1);
eta_node=coord(:,2);
zeta_node=coord(:,3);
Lmat=zeros(8,8);
%calculation of L transformation matrix
for count_i = 1:8
    for count_j = 1:size(P_g,1)
        Lmat(count_j,count_i) = 1/8*(1+xi_node(count_i)/P_g(count_j,1))*(1+eta_node(count_i)/P_g(count_j,2))*(1+zeta_node(count_i)/P_g(count_j,3));
    end
end


end % function end