function [KG,Melc,Vel,DPel,T,B,SG,Gel,Cel]=Penta(E,NU,rho,coord,n,Ns,Nt,Nu,GyroFlag)
% PENTA 6 ELEMENT FORMULATION
% decomposition of each penta in 3 hexa.
% Computes the required matrices for a 6 node Pentahedron splitting it
% in 3 hexahedra for the stiffness matrix, and using Ansys formulation
% for the Mass, Coriolis and Gyroscopic skew-symmetric matrices by
% considering the pentahedron as an hexadron with 2 degenerate nodes.
%
% INPUTS (E,NU,rho,coord)
% E: Material modulus of elasticity
% NU: Poisson ratio of the material
% rho: Material density
% coordor: Coordinates of the 6 nodes of the element
% GyroFlag: set to 1 if gyroscopic and coriolis matrices are needed, 0
% otherwise
% (Material properties are without dimensions since they are related to the
% metric used in the coordinates and therefore are user-dependent. For S.I.
% if the coordinates are in [mm], then E->[MPa], rho->[ton/mm^3]
% if the coordinates are in [m], then E->[Pa], rho->[kg/m^3])
%
% OUTPUTS[KG,Melc,Mell,Vel,DPel,btot,T,B,D,btotn,Gel,Cel]
% KG: Stiffness matrix (18x18)
% Melc: Consistent mass matrix (18x18)
% Mell: Lumped mass matrix (18x18 diagonal)
% Vel: Element volume
% DPel: Element distortion
% btot:
% T: Transformation matrix to take into account the enhanced assumed strain
% states when computing element strains and stresses
% B: Centroid Displacement/Strain relationship matrix
% D: Element material matrix
% btotn: Nodal Displacement/Strain relationship matrix
% Gel: Element Gyroscopic matrix (18x18 skew-symmetric) if requested
% Cel: Element Coriolis matrix (18x18 skew-symmetric) if requested


coord6=coord;

CenterA=mean(coord(1:3,:)); % centroid of first triangle face
CenterB=mean(coord(4:6,:)); % centroid of second triangle face
% additional needed points
% first face
PC=mean(coord(1:2,:));
PD=mean(coord(2:3,:));
PE=mean(coord([1,3],:));
% second face
PF=mean(coord(4:5,:));
PG=mean(coord(5:6,:));
PH=mean(coord([4,6],:));

% First Hexa: [1 C A E, 4 F B H]
% Second Hexa: [2 D A C, 5 G B F]
% Third Hexa: [3 E A D, 6 H B G]

elementi=[1 8 9 1 7 12 14 4 11;
    2 8 7 2 10 12 11 5 13;
    3 8 10 3 9 12 13 6 14];

doforder8=[1 9 17 2 10 18 3 11 19 4 12 20 5 13 21 6 14 22 7 15 23 8 16 24];

nodi=[1 coord(1,:);
    2 coord(2,:);
    3 coord(3,:);
    4 coord(4,:);
    5 coord(5,:);
    6 coord(6,:);
    7 PC;
    8 CenterA;
    9 PE;
    10 PD;
    11 PF;
    12 CenterB;
    13 PG;
    14 PH];


nnodi=length(nodi);

% creazione mappa elementi con numerazione nodi matlab
IDnodi=nodi(:,1);
nele=3;
elemconv=zeros(nele,8);
Kass=zeros(nnodi*3,nnodi*3);
Ktemp=zeros(24,24,3);
Kel=zeros(24,24);

% nella numerazione Nastran i nodi non partono da 1 e segue una
    % logica diversa, con questa linea a ogni nodo si associa la riga
    % della matrice matlab in cui sono presenti le coordinate
    [~,elemconv(:,:)] = builtin('_ismemberhelper',elementi(:,2:9),IDnodi(:,1));
    
coord8=zeros(8,3);
% Stiffness computation
for ii=1:3

    % genera la matrice delle coordinate locali pescando dalla matrice
    % delle coordinate dei nodi
    coord8(1:8,1:3)=nodi(elemconv(ii,:),2:4);
    
    % calcolo matrici dell'elemento
    [Kel,~,~,~,~,~,~,~,~,Delm,~]=Hexa(E,NU,rho,coord8,n,Ns,Nt,Nu,GyroFlag);
    Ktemp(:,:,ii)=Kel;
    % i gdl totali sono 3* numero nodi
    % nel sistema di riferimento locale i gdl sono
    % [x1,x2,....,y1,y2,....,z1,z2,.....] mentre nel sistema globale
    % [x1,y1,z1,....xn,yn,zn,....]
    %mappasl=[elemconv(ii,:)*3-2,elemconv(ii,:)*3-1,elemconv(ii,:)*3];
    %%
    mappasl=repmat(elemconv(ii,:),1,3);
    mappasl=mappasl(1,doforder8)*3+repmat([-2:1:0],1,8);
%     map(ii,:)=mappasl;
    Kass(mappasl,mappasl)=Kass(mappasl,mappasl)+Kel(1:24,1:24);
    %
end

%% HEXA assembly
% row=reshape(repmat(map,1,24)',24,24,3);
% col=permute(row,[2 1 3]);
% Kass=sparse(col(:),row(:),Ktemp(:));

% MASTER - SLAVE MPC
lat=1/2;
cent=1/3;
rbar=zeros(42,18);
rbar(1:18,1:18)=eye(18);
% center A N8
rbar(22:24,1:9)=cent*[eye(3),eye(3),eye(3)];
% center B N12
rbar(34:36,10:18)=cent*[eye(3),eye(3),eye(3)];
% pC
rbar(19:21,1:6)=lat*[eye(3),eye(3)];
% pD
rbar(28:30,4:9)=lat*[eye(3),eye(3)];
% pE
rbar(25:27,[1:3,7:9])=lat*[eye(3),eye(3)];
% pF
rbar(31:33,10:15)=lat*[eye(3),eye(3)];
% pG
rbar(37:39,13:18)=lat*[eye(3),eye(3)];
% pH
rbar(40:42,[10:12,16:18])=lat*[eye(3),eye(3)];
% MASTER - SLAVE reduction
KG=rbar'*Kass*rbar;

% Inertial matrices computation
% degenerated hexa
coordor=[coord6(1:3,:);coord6(3,:);coord6(4:6,:);coord6(6,:)];
[MC8,Vel,DPel,Ge8,Ce8]=Masses(rho,coordor,n,Ns,Nt,Nu,GyroFlag);

MC8=reorder_dof(MC8);
% condense the values of the degenerated nodes to the existing ones
[Melc]=condense(MC8);

if GyroFlag == 1
    [Gel]=condense(Ge8);
    [Cel]=condense(Ce8);
    Gel=Gel-Gel.';
else
    Gel=0;
    Cel=0;
end

% displacement- strain/stress matrices averaging and composition
[btotg,detJt,~,detrap, B]=Bmatrix(coordor,Ns,Nt,Nu);
[~,T]= Kbrick(btotg,Delm,detJt,detrap);
SG=zeros(48,30);
SG(1:6,:)=Delm*btotg(:,:,1);
SG(7:12,:)=Delm*btotg(:,:,2);
SG(13:18,:)=Delm*btotg(:,:,3);
SG(19:24,:)=Delm*btotg(:,:,4);
SG(25:30,:)=Delm*btotg(:,:,5);
SG(31:36,:)=Delm*btotg(:,:,6);
SG(37:42,:)=Delm*btotg(:,:,7);
SG(43:48,:)=Delm*btotg(:,:,8);

B=B/Vel;

end %function end
