function [n27,Ns27,Nt27,Nu27,Ns8,Nt8,Nu8]= ShapeFunsH20

n27=zeros(3,60,27);
Ns27=zeros(27,20);
Nt27=zeros(27,20);
Nu27=zeros(27,20);

% 3x3x3
GaussP=[-0.774596669241483,-0.774596669241483,-0.774596669241483;
    0,-0.774596669241483,-0.774596669241483;
    0.774596669241483,-0.774596669241483,-0.774596669241483;
    -0.774596669241483,0,-0.774596669241483;
    0,0,-0.774596669241483;
    0.774596669241483,0,-0.774596669241483;
    -0.774596669241483,0.774596669241483,-0.774596669241483;
    0,0.774596669241483,-0.774596669241483;
    0.774596669241483,0.774596669241483,-0.774596669241483;
    -0.774596669241483,-0.774596669241483,0;
    0,-0.774596669241483,0;
    0.774596669241483,-0.774596669241483,0;
    -0.774596669241483,0,0;
    0,0,0;
    0.774596669241483,0,0;
    -0.774596669241483,0.774596669241483,0;
    0,0.774596669241483,0;
    0.774596669241483,0.774596669241483,0;
    -0.774596669241483,-0.774596669241483,0.774596669241483;
    0,-0.774596669241483,0.774596669241483;
    0.774596669241483,-0.774596669241483,0.774596669241483;
    -0.774596669241483,0,0.774596669241483;
    0,0,0.774596669241483;
    0.774596669241483,0,0.774596669241483;
    -0.774596669241483,0.774596669241483,0.774596669241483;
    0,0.774596669241483,0.774596669241483;
    0.774596669241483,0.774596669241483,0.774596669241483];

for ii=1:27
    
      
    xi=GaussP(ii,1);
    et=GaussP(ii,2);
    ze=GaussP(ii,3);
    
    omg=1-xi;
    omh=1-et;
    omr=1-ze;
    opg=1+xi;
    oph=1+et;
    opr=1+ze;
    tpgphpr=opg+oph+ze;
    tmgphpr=omg+oph+ze;
    tmgmhpr=omg+omh+ze;
    tpgmhpr=opg+omh+ze;
    tpgphmr=opg+oph-ze;
    tmgphmr=omg+oph-ze;
    tmgmhmr=omg+omh-ze;
    tpgmhmr=opg+omh-ze;
    omgopg=omg*opg/4;
    omhoph=omh*oph/4;
    omropr=omr*opr/4;
    omgmopg=(omg-opg)/4;
    omhmoph=(omh-oph)/4;
    omrmopr=(omr-opr)/4;
    
    %     shape functions
    
    N(1,1)=-omg*omh*omr*tpgphpr/8;
    N(1,2)=-opg*omh*omr*tmgphpr/8;
    N(1,3)=-opg*oph*omr*tmgmhpr/8;
    N(1,4)=-omg*oph*omr*tpgmhpr/8;
    N(1,5)=-omg*omh*opr*tpgphmr/8;
    N(1,6)=-opg*omh*opr*tmgphmr/8;
    N(1,7)=-opg*oph*opr*tmgmhmr/8;
    N(1,8)=-omg*oph*opr*tpgmhmr/8;
    N(1,9)=omgopg*omh*omr;
    N(1,10)=omhoph*opg*omr;
    N(1,11)=omgopg*oph*omr;
    N(1,12)=omhoph*omg*omr;
    N(1,13)=omgopg*omh*opr;
    N(1,14)=omhoph*opg*opr;
    N(1,15)=omgopg*oph*opr;
    N(1,16)=omhoph*omg*opr;
    N(1,17)=omropr*omg*omh;
    N(1,18)=omropr*opg*omh;
    N(1,19)=omropr*opg*oph;
    N(1,20)=omropr*omg*oph;
    
    % matrice funzioni di forma
    n27(:,:,ii)=[N,zeros(1,40)
        zeros(1,20),N,zeros(1,20)
        zeros(1,40),N];
    
    % local derivatives of the shape functions: xi-derivative
    
    Nxi(1, 1)=omh*omr*(tpgphpr-omg)/8;
    Nxi(1, 2)=(opg-tmgphpr)*omh*omr/8;
    Nxi(1, 3)=(opg-tmgmhpr)*oph*omr/8;
    Nxi(1, 4)=oph*omr*(tpgmhpr-omg)/8;
    Nxi(1, 5)=omh*opr*(tpgphmr-omg)/8;
    Nxi(1, 6)=(opg-tmgphmr)*omh*opr/8;
    Nxi(1, 7)=(opg-tmgmhmr)*oph*opr/8;
    Nxi(1, 8)=oph*opr*(tpgmhmr-omg)/8;
    Nxi(1, 9)=omgmopg*omh*omr;
    Nxi(1,10)=omhoph*omr;
    Nxi(1,11)=omgmopg*oph*omr;
    Nxi(1,12)=-omhoph*omr;
    Nxi(1,13)=omgmopg*omh*opr;
    Nxi(1,14)=omhoph*opr;
    Nxi(1,15)=omgmopg*oph*opr;
    Nxi(1,16)=-omhoph*opr;
    Nxi(1,17)=-omropr*omh;
    Nxi(1,18)=omropr*omh;
    Nxi(1,19)=omropr*oph;
    Nxi(1,20)=-omropr*oph;
    
    %   local derivatives of the shape functions: eta-derivative
    
    Net(1, 1)=omg*omr*(tpgphpr-omh)/8;
    Net(1, 2)=opg*omr*(tmgphpr-omh)/8;
    Net(1, 3)=opg*(oph-tmgmhpr)*omr/8;
    Net(1, 4)=omg*(oph-tpgmhpr)*omr/8;
    Net(1, 5)=omg*opr*(tpgphmr-omh)/8;
    Net(1, 6)=opg*opr*(tmgphmr-omh)/8;
    Net(1, 7)=opg*(oph-tmgmhmr)*opr/8;
    Net(1, 8)=omg*(oph-tpgmhmr)*opr/8;
    Net(1, 9)=-omgopg*omr;
    Net(1,10)=omhmoph*opg*omr;
    Net(1,11)=omgopg*omr;
    Net(1,12)=omhmoph*omg*omr;
    Net(1,13)=-omgopg*opr;
    Net(1,14)=omhmoph*opg*opr;
    Net(1,15)=omgopg*opr;
    Net(1,16)=omhmoph*omg*opr;
    Net(1,17)=-omropr*omg;
    Net(1,18)=-omropr*opg;
    Net(1,19)=omropr*opg;
    Net(1,20)=omropr*omg;
    
    %     local derivatives of the shape functions: zeta-derivative
    
    Nze(1, 1)=omg*omh*(tpgphpr-omr)/8;
    Nze(1, 2)=opg*omh*(tmgphpr-omr)/8;
    Nze(1, 3)=opg*oph*(tmgmhpr-omr)/8;
    Nze(1, 4)=omg*oph*(tpgmhpr-omr)/8;
    Nze(1, 5)=omg*omh*(opr-tpgphmr)/8;
    Nze(1, 6)=opg*omh*(opr-tmgphmr)/8;
    Nze(1, 7)=opg*oph*(opr-tmgmhmr)/8;
    Nze(1, 8)=omg*oph*(opr-tpgmhmr)/8;
    Nze(1, 9)=-omgopg*omh;
    Nze(1,10)=-omhoph*opg;
    Nze(1,11)=-omgopg*oph;
    Nze(1,12)=-omhoph*omg;
    Nze(1,13)=omgopg*omh;
    Nze(1,14)=omhoph*opg;
    Nze(1,15)=omgopg*oph;
    Nze(1,16)=omhoph*omg;
    Nze(1,17)=omrmopr*omg*omh;
    Nze(1,18)=omrmopr*opg*omh;
    Nze(1,19)=omrmopr*opg*oph;
    Nze(1,20)=omrmopr*omg*oph;
    
    Ns27(ii,:)=Nxi;
    Nt27(ii,:)=Net;
    Nu27(ii,:)=Nze;
    
end

% 2x2x2
Ns8=zeros(8,20);
Nt8=zeros(8,20);
Nu8=zeros(8,20);

pos=sqrt(3)/3;

GaussP=[-pos,-pos,-pos;
        pos,-pos,-pos;
        pos,pos,-pos;
        -pos,pos,-pos;
        -pos,-pos,pos;
        pos,-pos,pos;
        pos,pos,pos;
        -pos,pos,pos];
    
for ii=1:8
    
    xi=GaussP(ii,1);
    et=GaussP(ii,2);
    ze=GaussP(ii,3);
    
    omg=1-xi;
    omh=1-et;
    omr=1-ze;
    opg=1+xi;
    oph=1+et;
    opr=1+ze;
    tpgphpr=opg+oph+ze;
    tmgphpr=omg+oph+ze;
    tmgmhpr=omg+omh+ze;
    tpgmhpr=opg+omh+ze;
    tpgphmr=opg+oph-ze;
    tmgphmr=omg+oph-ze;
    tmgmhmr=omg+omh-ze;
    tpgmhmr=opg+omh-ze;
    omgopg=omg*opg/4;
    omhoph=omh*oph/4;
    omropr=omr*opr/4;
    omgmopg=(omg-opg)/4;
    omhmoph=(omh-oph)/4;
    omrmopr=(omr-opr)/4;
    
    % local derivatives of the shape functions: xi-derivative
    
    Nxi(1, 1)=omh*omr*(tpgphpr-omg)/8;
    Nxi(1, 2)=(opg-tmgphpr)*omh*omr/8;
    Nxi(1, 3)=(opg-tmgmhpr)*oph*omr/8;
    Nxi(1, 4)=oph*omr*(tpgmhpr-omg)/8;
    Nxi(1, 5)=omh*opr*(tpgphmr-omg)/8;
    Nxi(1, 6)=(opg-tmgphmr)*omh*opr/8;
    Nxi(1, 7)=(opg-tmgmhmr)*oph*opr/8;
    Nxi(1, 8)=oph*opr*(tpgmhmr-omg)/8;
    Nxi(1, 9)=omgmopg*omh*omr;
    Nxi(1,10)=omhoph*omr;
    Nxi(1,11)=omgmopg*oph*omr;
    Nxi(1,12)=-omhoph*omr;
    Nxi(1,13)=omgmopg*omh*opr;
    Nxi(1,14)=omhoph*opr;
    Nxi(1,15)=omgmopg*oph*opr;
    Nxi(1,16)=-omhoph*opr;
    Nxi(1,17)=-omropr*omh;
    Nxi(1,18)=omropr*omh;
    Nxi(1,19)=omropr*oph;
    Nxi(1,20)=-omropr*oph;
    
    %   local derivatives of the shape functions: eta-derivative
    
    Net(1, 1)=omg*omr*(tpgphpr-omh)/8;
    Net(1, 2)=opg*omr*(tmgphpr-omh)/8;
    Net(1, 3)=opg*(oph-tmgmhpr)*omr/8;
    Net(1, 4)=omg*(oph-tpgmhpr)*omr/8;
    Net(1, 5)=omg*opr*(tpgphmr-omh)/8;
    Net(1, 6)=opg*opr*(tmgphmr-omh)/8;
    Net(1, 7)=opg*(oph-tmgmhmr)*opr/8;
    Net(1, 8)=omg*(oph-tpgmhmr)*opr/8;
    Net(1, 9)=-omgopg*omr;
    Net(1,10)=omhmoph*opg*omr;
    Net(1,11)=omgopg*omr;
    Net(1,12)=omhmoph*omg*omr;
    Net(1,13)=-omgopg*opr;
    Net(1,14)=omhmoph*opg*opr;
    Net(1,15)=omgopg*opr;
    Net(1,16)=omhmoph*omg*opr;
    Net(1,17)=-omropr*omg;
    Net(1,18)=-omropr*opg;
    Net(1,19)=omropr*opg;
    Net(1,20)=omropr*omg;
    
    %     local derivatives of the shape functions: zeta-derivative
    
    Nze(1, 1)=omg*omh*(tpgphpr-omr)/8;
    Nze(1, 2)=opg*omh*(tmgphpr-omr)/8;
    Nze(1, 3)=opg*oph*(tmgmhpr-omr)/8;
    Nze(1, 4)=omg*oph*(tpgmhpr-omr)/8;
    Nze(1, 5)=omg*omh*(opr-tpgphmr)/8;
    Nze(1, 6)=opg*omh*(opr-tmgphmr)/8;
    Nze(1, 7)=opg*oph*(opr-tmgmhmr)/8;
    Nze(1, 8)=omg*oph*(opr-tpgmhmr)/8;
    Nze(1, 9)=-omgopg*omh;
    Nze(1,10)=-omhoph*opg;
    Nze(1,11)=-omgopg*oph;
    Nze(1,12)=-omhoph*omg;
    Nze(1,13)=omgopg*omh;
    Nze(1,14)=omhoph*opg;
    Nze(1,15)=omgopg*oph;
    Nze(1,16)=omhoph*omg;
    Nze(1,17)=omrmopr*omg*omh;
    Nze(1,18)=omrmopr*opg*omh;
    Nze(1,19)=omrmopr*opg*oph;
    Nze(1,20)=omrmopr*omg*oph;
    
    Ns8(ii,:)=Nxi;
    Nt8(ii,:)=Net;
    Nu8(ii,:)=Nze;
    
end

end %function end