function [Kass,Mass,Massl,Gyro,Corio,vtot,distorsione,Te,BK,BGnew,SG,map,elemconv,nodi2,elementi2,ListPenta]=poliFEMo_sparse(nodi,elementi,E,NU,rho,GyroFlag,order,elem_material)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%                               poliFEMo
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Program to compute the structural matrices from a finite element mesh.
%
% INPUT DESCRIPTION
%
% REQUIRED INPUTS
%[outputs]=poliFEMo_sparse(nodi,elementi,E,NU,rho)
% ALL INPUTS
%[outputs]=poliFEMo_sparse(nodi,elementi,E,NU,rho,GyroFlag,order,elem_material)
%
% -------------------------------------------------------------------------
% nodi = [NodeNumber Xcoo Ycoo Zcoo] - size(NumberOfNodes,4)
% Finite element mesh nodes number and coordinates in a global cartesian
% reference frame
%
% Caution: The nodes MUST be sorted in increasing order. If not the program
% will throw a WARNING and will automatically sort the nodes in increasing
% order. In this case to correctly use the map and elemconv outputs the
% sorted list must be used and the program provides this list as outputs
% nodi2 and elementi2.
%
% -------------------------------------------------------------------------
% elementi = [ElementNumber Node1 Node2 ... NodeN] - size(NumberOfElements,9 or 21)
% Finite elements of the mesh. The node numbering follows the convention of
% Ansys for both first and second order elements. 6 or 15 nodes wedge
% (Penta6 or Penta15) elements must have repeated indices. For first order
% elements the element matrix must have 9 colums in total, while it must
% have 21 columns for second order elements. For the node numbering
% convention please refer to the following links for both hexa and penta
% elements.
% First order:
% https://www.sharcnet.ca/Software/Ansys/16.2.3/en-us/help/ans_elem/Hlp_E_SOLID185.html
% [ElemNumber I J K L M N O P] Hexa8
% [ElemNumber I J K K M N O O] Penta6
% Second order:
% https://www.sharcnet.ca/Software/Ansys/16.2.3/en-us/help/ans_elem/Hlp_E_SOLID186.html
% [ElemNumber I J K L M N O P Q R S T U V W X Y Z A B] Hexa20
% [ElemNumber I J K K M N O O Q R K T U V O X Y Z A A] Penta15
%
% -------------------------------------------------------------------------
% E = [Young's modulus] -  size(NumberOfMaterials,1)
% Material modulus of elasticity. More than one materials can be specified
% for the model, but one element can have only one material, so the number
% of materials can range from 1 to NumberOfElements. NumberOfMaterials must
% be the same for each material property (E,NU,rho).
%
% -------------------------------------------------------------------------
% NU = [Poisson Ratio] -  size(NumberOfMaterials,1)
% Material Poisson ratio. More than one materials can be specified
% for the model, but one element can have only one material, so the number
% of materials can range from 1 to NumberOfElements. NumberOfMaterials must
% be the same for each material property (E,NU,rho).
%
% -------------------------------------------------------------------------
% rho = [density] -  size(NumberOfMaterials,1)
% Material density. More than one materials can be specified
% for the model, but one element can have only one material, so the number
% of materials can range from 1 to NumberOfElements. NumberOfMaterials must
% be the same for each material property (E,NU,rho).
%
% -------------------------------------------------------------------------
% GyroFlag= [flag] - scalar( 0 or 1 ) - Default: 0 - OPTIONAL
% Flag to request or not the computation of the Gyroscopic and Coriolis
% matrices to include rotational effects (around the global Z axis)
% in a modal analysis. Input 1 to request the matrices
%  or leave the field to 0 or blank otherwise. For an
% analysis in the rotating frame using the Coriolis matrix, remember to
% include the so called "spin softening" effect by issuing:
% SpinSoftening= MassMatrix;
% SpinSoftening(3:3:end,3:3:end)=0;
%
% -------------------------------------------------------------------------
% order = [order of the elements] - scalar( 1 or 2 ) - Default: Autodefined
% - OPTIONAL
% Flag to specify the order of the elements to be used, optionally set it
% to 1 to force first order elements or 2 to force second order elements.
% If not specified the program automatically determines the order of the
% elements by the size of the element matrix (9 for first order and 21 for
% second order). Not yet implemented: If specified 1 and the element matrix
% is of the second order the program will automatically "reduce" the
% elements to first order by ignoring the midside nodes. The map and
% elemeconv outputs are redimensioned accordingly and the nodi2 and
% elementi2 will reflect this "reduction" of order. On the opposite, if
% specified 2 and the element matrix is of the first order the program will
% automatically remesh the model to add the midside nodes. The map and
% elemeconv outputs are redimensioned accordingly and the nodi2 and
% elementi2 will reflect this remeshing of finite element mesh.
%
% -------------------------------------------------------------------------
% elem_material = [ElementNumber MaterialNumber] - size(NumberOfElements,2)
% - OPTIONAL (REQUIRED if NumberOfMaterials is greater than 1)
% Table to associate each element with its own material. 
% The program for element I, will look in the elem_material correspondig
% to element I the MaterialNumber K and will use the values of E,NU and rho
% listed in row K and will use these values to compute the structural
% matrices required.
%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUT DESCRIPTION
% [Kass,Mass,Massl,Gyro,Corio,vtot,distorsione,Te,BK,BG,SG,map,elemconv,nodi2,elementi2,ListPenta]
%
%
% -------------------------------------------------------------------------
% Kass = [Stiffness matrix] - size(NumberOfNodes*3,NumberOfNodes*3)
% Assembled structural stiffness matrix
%
% -------------------------------------------------------------------------
% Mass = [Consistent mass matrix] - size(NumberOfNodes*3,NumberOfNodes*3)
% Assembled consistent mass matrix
%
% -------------------------------------------------------------------------
% Massl = [Lumped mass matrix] - size(NumberOfNodes*3,NumberOfNodes*3)
% Assembled lumped mass matrix
%
% -------------------------------------------------------------------------
% Gyro = [Gyroscopic matrix] - size(NumberOfNodes*3,NumberOfNodes*3)
% Assembled gyroscopic effect matrix. Used for analysis of components
% rotating around the global Z axis in a Stationary reference frame.
% Optional matrix: issued only if GyroFlag = 1, otherwise the output is
% 'Not Requested'
%
% -------------------------------------------------------------------------
% Corio = [Coriolis matrix] - size(NumberOfNodes*3,NumberOfNodes*3)
% Assembled Coriolis effect matrix. Used for analysis of components
% rotating around the global Z axis in a Rotating reference frame. For an
% analysis in the rotating frame using the Coriolis matrix, remember to
% include the so called "spin softening" effect by issuing:
% SpinSoftening= MassMatrix;
% SpinSoftening(3:3:end,3:3:end)=0;
% Optional matrix: issued only if GyroFlag = 1, otherwise the output is
% 'Not Requested'
%
% -------------------------------------------------------------------------
% vtot = [model volume] - size(1,2)
% Volume of the model. The first column value indicates the total volume of
% the model regardless of the element type, while the second column
% indicates the volume occupied by wedge (Penta6 or Penta15) elements.
%
% -------------------------------------------------------------------------
% distorsione = [Element distortion index] - size(NumberOfElements,2)
% Distortion index of each element. Lower values indicate low distorstion
% (good). Too distorted elements can cause unrealistic results. Useful for
% mesh quality analysis.
%
% -------------------------------------------------------------------------
% Te = [Transformation matrix] - size(6,24,NumberOfElements) (or empty)
% Matrices for each element to compute the additional hidden displacements
% of Hexa8 elements. Only issued if elements of this type are present in
% the model, otherwise it's empty. Used in stress recovery.
%
% -------------------------------------------------------------------------
% BK = [Strain displacement at centroid] - size(variable depending on
% element type,variable depending on element type, NumberOfElements)
% Strain-displacement matrix evaluated at element centroid, used for strain
% recovery at centroid location of each element.
%
% -------------------------------------------------------------------------
% BG = [Strain-displacement at integration] - size(variable depending on
% element type,variable depending on element type, NumberOfElements)
% Strain-displacement matrix evaluated at element integration points 
% location (Gauss points), used for strain recovery at centroid location of
% each element.
%
% -------------------------------------------------------------------------
% SG = [Stress-displacement at integration] - size(variable depending on
% element type,variable depending on element type, NumberOfElements)
% Stress-displacement matrix evaluated at element integration points 
% location (Gauss points), used for stress recovery at centroid location of
% each element.
%
% -------------------------------------------------------------------------
% map = [dof map] - size([24 or 60], NumberOfElements)
% Degrees of freedom location of each node in the structural matrices, used
% in conjunction with output elemconv to correctly pinpoint the degrees of
% freedom associated to a certain node inside the structural matrices. The 
% row size is 24 for first order elements and 60 for the second one. Penta
% (6 or 15) elements will have zeroes in their row since the duplicated
% values are removed.
%
% -------------------------------------------------------------------------
% elemconv = [node row map] - size([8 or 20], NumberOfElements)
% Node rows location of each node in the structural matrices, used
% in conjunction with output map to correctly pinpoint the degrees of
% freedom associated to a certain node inside the structural matrices. Each
% value in this matrix indicated the row associated to the corresponding
% node of the element. The row size is 8 for first order elements or 20 for
% the second one. Penta (6 or 15) will have zeroes in their row since the
% duplicated values are removed.
%
% To correctly find the dofs associated to a certain node, one could use
% for example:
%    % Find target node NodeTarg
%      [rn,~]=find(NodeTarg==nodi(:,1));
%      [row,col]=find(rn(1)==elemconv);
%      pos=map(row(1),3*(col(1)-1)+1:3*(col(1)-1)+3); %XYZ dof
%    % the values in  pos indicate the position of the X Y and Z dofs
%    % associated to node NodeTarg
%
% -------------------------------------------------------------------------
% nodi2 = [sorted or modified nodes] - size(variable,4)
% Sorted node list or node list of the modified nodes. If the input node
% list (nodi) was not input by the user in sorted order the program will
% automatically sort them and output the sorted list in this variable. This
% sorted list will have to be used in conjuction to map and elemconv to
% correctly find the dofs associated to a certain node. The program issues
% a warning in this case. Insted if the user inputs a mesh not
% corresponding to the order of the analysis selected (optional input
% order) the program will automatically remesh the model to correspond to
% the selected order. The node list in this output will reflect these
% changes in this case. For example if a second order mesh (with midside 
% nodes) was input but order =1 , then the program will automatically
% remove the additional midside nodes and the node list in this output will
% contain only the corner nodes. Viceversa if a first order mesh was input
% but a second order analysis was selected (order = 2) the the program will
% add all required midside nodes and this output will contain the complete
% list of the nodes used.
%
% -------------------------------------------------------------------------
% elementi2 = [modified elements] - size(NumberOfElements,9 or 21)
% Modified element list and connectivity. If the user inputs a mesh not
% corresponding to the order of the analysis selected (optional input
% order) the program will automatically remesh the model to correspond to
% the selected order. The element list in this output will reflect these
% changes in this case. For example if a second order mesh (with midside 
% nodes) was input but order =1 , then the program will automatically
% remove the additional midside nodes and the element list in this output 
% will contain only the corner nodes. Viceversa if a first order mesh was
% input but a second order analysis was selected (order = 2) the the 
% program will add all required midside nodes and this output will contain 
% the complete list of the elements used.
%
% -------------------------------------------------------------------------
% ListPenta = [list of Penta elements] - size(variable,2)
% List of Penta (6 or 15) contained in the finite element model. Each row
% (if any) contains in the first column the element number corresponding to
% the wedge element as present in the elementi input table. In the second
% column is indicated the row correspondig to said element for easier
% identification in case of oddly numbered elements. If this output is
% equal to [0 0] it means that no wedge elements are present in the input
% mesh.
%
% -------------------------------------------------------------------------
% BETA v0.9 F. Bruzzone 12/2018


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INIT
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% conversione numeri nodi
nnodi=length(nodi);
nele=size(elementi,1);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CHECKS
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% check iniziali su numero variabili in input, minimo  5
% nodi,elementi,E,NU,rho
narginchk(5,8);

if nargin <= 5
    multimaterial=0;
    GyroFlag =0;
    order= 1;
    nodi2=0;
    elementi2=0;
end
if nargin <= 6
    multimaterial=0;
    order = 1;
    nodi2=0;
    elementi2=0;
end
if nargin <= 7
    multimaterial=0;
end

if size(E,1) > 2 ||  size(NU,1) > 2 || size(rho,1) > 2
    multimaterial=1;
else
    multimaterial = 0;
    Ei=E(1);
    NUi=NU(1);
    rhoi=rho(1);
end

% node numbering check
if ~issorted(nodi(:,1),1,'ascend')
    if nargout<=13
        warning('WARNING: The node input supplied is not in sorted order. The program will automatically sort them in the correct order, but the sorted node list will not be output since it was not required. Restart the analysis by first sorting the nodes in ascending order, or enable the last two outputs in order to be able to correctly use the element conversion table and the d.o.f.s map. Press F5 or click continue to resume the analysis or abort.')
        keyboard
    else
        warning('WARNING: The node input supplied is not in sorted order. The program will automatically sort them in the correct order, remember to use the sorted (output 14) node list to correctly use the element conversion table and the d.o.f.s map.')
    end
    
    [~,idx]=sort(nodi(:,1),1,'ascend');
    ns=nodi(idx,:);
    clear nodi
    nodi=ns;
    
end

% elem_material= [N° elemento, N° materiale] size-> (nele,2)
% check che le dimensioni dei vettori del materiale siano corrette
if multimaterial == 1
    if size(elem_material,1) ~= nele
        error('ERROR: More than 1 material data were specified as input, but the size of the vector containing the Elem N° <--> Material N° is not of the correct size. Check the input data.')
    end
    
    if max(elem_material(:,2)) > size(E,1)
        error('ERROR: More than 1 material data were specified as input, but one of the elements points to a E value that is out of range. Check that the number of E values specified corresponds to the number of materials specified for the elements.')
    end
    
    if max(elem_material(:,2)) > size(NU,1)
        error('ERROR: More than 1 material data were specified as input, but one of the elements points to a NU value that is out of range. Check that the number of NU values specified corresponds to the number of materials specified for the elements.')
    end
    if max(elem_material(:,2)) > size(rho,1)
        error('ERROR: More than 1 material data were specified as input, but one of the elements points to a rho value that is out of range. Check that the number of rho values specified corresponds to the number of materials specified for the elements.')
    end
    
    if multimaterial == 1 && ~exist('elem_material','var')
        warning('WARNING: more than one material data were specified, but the element-material correspondence were not. If you choose to continue all elements will use the first material.');
        keyboard
        multimaterial = 0;
        E=Ei(1);
        NU=NUi(1);
        rho=rhoi(1);
    end
    
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% nel caso sia richiesto, uso elementi del secondo ordine
if order == 2 && size(elementi,2)<10
    error('buggato, stop'); % per sicurezza
    keyboard 
    [nodi2,elementi2]=RemeshTo2ndOrderElements(nodi,elementi);
elseif order == 2 && size(elementi,2)==21
    nodi2=nodi;
    elementi2=elementi;
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FIRST ORDER ELEMENTS ----------------------------------------------------
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% creazione mappa elementi con numerazione nodi matlab
if order == 1
    % dummy per output
    nodi2=nodi;
    elementi2=elementi;
    
    % calcolo funzioni di forma e relative derivate ai punti di
    % integrazione
    [n,Ns,Nt,Nu]=ShapeFuns8;
    
    IDnodi=nodi(:,1);
    elemconv=zeros(nele,8);
    distorsione=zeros(nele,2);
    numhexa=0;
    numpenta=0;
    
    if GyroFlag==0
        Gyro= 'Not Requested';
        Corio='Not Requested';
    end
    
    % preallocate matrices
    % Strain displacement matrix at element centroid
    BK=zeros(6,30,nele);
    % Strain displacement matrices at gauss point locations
    BG=zeros(48,30,nele);
    % Strain displacement matrices at gauss point locations -
    % Newformulation
    BGnew=zeros(6,24,8,nele);
    % Stress displacement matrices at gauss point locations
    SG=zeros(48,30,nele);
    map=zeros(nele,24);
    
    % transformation matrices
    Te=zeros(6,24,nele);
    
    % Volume
    % vtot = [ Total Volume (hexa+penta), Penta only volume];
    vtot=[0 0];
    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];
    doforder6=[1:6:18 2:6:18 3:6:18 4:6:18 5:6:18 6:6:18];
    
    % temporary matrices
    [Hexanum,Pentanum,ListPenta]=GetElemTypes(elementi);
    % penta
    KtempP=zeros(18,18,Pentanum);
    MtempP=zeros(18,18,Pentanum);
    MLtempP=zeros(18,18,Pentanum);
    
    if GyroFlag==1
        GtempP=zeros(18,18,Pentanum);
        CtempP=zeros(18,18,Pentanum);
    end
    map_penta=zeros(Pentanum,18);
    listpenta=zeros(Pentanum,1);
    
    %hexa
    Ktemp=zeros(24,24,Hexanum);
    Mtemp=zeros(24,24,Hexanum);
    MLtemp=zeros(24,24,Hexanum);
    if GyroFlag==1
        Gtemp=zeros(24,24,Hexanum);
        Ctemp=zeros(24,24,Hexanum);
    end
    map_hexa=zeros(Hexanum,24);
    listhexa=zeros(Hexanum,1);
    
%     [~,elemconv(:,:)] = builtin('_ismemberhelper',elementi(:,2:9),IDnodi(:,1));
        [~,elemconv(:,:)] = ismember(elementi(:,2:9),IDnodi(:,1));
    
    % Element assembly loop
    for ii=1:nele
      
        % genera la matrice delle coordinate locali pescando dalla matrice
        % delle coordinate dei nodi
        coord(1:8,1:3)=nodi(elemconv(ii,:),2:4);
        
        pentaflag=0; % 0 Hexa 1 Penta
        
        if elemconv(ii,3) == elemconv(ii,4) && elemconv(ii,7) == elemconv(ii,8)
            % penta element
            pentaflag=1; % is penta
            
            % adjust elemconv table
            elemconv(ii,1:8)=[elemconv(ii,[1:3,5:7]), 0,0];
            
            % adjust coordinates
            coord6=coord([1:3,5:7],:);
            
            if multimaterial == 1
                Ei=E(elem_material(ii,2));
                NUi=NU(elem_material(ii,2));
                rhoi=rho(elem_material(ii,2));
            end
            
            [Kel,Melc,Vel,DPel,T,B,Stotg,Gel,Cel]=Penta(Ei,NUi,rhoi,coord6,n,Ns,Nt,Nu,GyroFlag);
            numpenta=numpenta+1;
            KtempP(:,:,numpenta)=Kel;
            MtempP(:,:,numpenta)=Melc;
            
            if GyroFlag==1
                GtempP(:,:,numpenta)=Gel;
                CtempP(:,:,numpenta)=Cel;
            end
            
            BK(:,:,ii)=B;
            SG(:,:,ii)=Stotg;
            Te(:,:,ii)=T;
            vtot=vtot+[Vel, Vel];
            distorsione(ii,:)=[elementi(ii,1),DPel];
            
          
        elseif elemconv(ii,3) ~= elemconv(ii,4) && elemconv(ii,7) ~= elemconv(ii,8)
            if multimaterial == 1
                Ei=E(elem_material(ii,2));
                NUi=NU(elem_material(ii,2));
                rhoi=rho(elem_material(ii,2));
            end
            
            % calcolo matrici dell'elemento
            [Kel,Melc,Gel,Cel,Vel,DPel,btotg,T,B,Delm,~,Bn]=Hexa(Ei,NUi,rhoi,coord,n,Ns,Nt,Nu,GyroFlag);
            numhexa=numhexa+1;
            Ktemp(:,:,numhexa)=Kel;
            Mtemp(:,:,numhexa)=Melc;
            if GyroFlag==1
                Gtemp(:,:,numhexa)=Gel;
                Ctemp(:,:,numhexa)=Cel;
            end
            
            BK(:,:,ii)=B;
            % New strain displacement calculation method
            BGnew(:,:,:,ii)=Bn;
            SG(1:6,:,ii)=Delm*btotg(:,:,1);
            SG(7:12,:,ii)=Delm*btotg(:,:,2);
            SG(13:18,:,ii)=Delm*btotg(:,:,3);
            SG(19:24,:,ii)=Delm*btotg(:,:,4);
            SG(25:30,:,ii)=Delm*btotg(:,:,5);
            SG(31:36,:,ii)=Delm*btotg(:,:,6);
            SG(37:42,:,ii)=Delm*btotg(:,:,7);
            SG(43:48,:,ii)=Delm*btotg(:,:,8);
            Te(:,:,ii)=T;
            vtot(1)=vtot(1)+Vel;
            distorsione(ii,:)=[elementi(ii,1),DPel];
        else
            Numerr=num2str(ii);
            error(['Element',' ', Numerr,' contains a number of nodes that does ', ...
                ' ','not correspond to an Hexahedron element (8 nodes) or a Pentahedron', ...
                ' element (6 nodes). Please review your mesh.']);
            
        end
        
        if pentaflag==0
            %% Hexa 
            % 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=repmat(elemconv(ii,:),1,3);
            mappasl=mappasl(1,doforder8)*3+repmat([-2:1:0],1,8);
            map(ii,:)=mappasl;
            
            map_hexa(numhexa,:)=mappasl;
            listhexa(numhexa,1)=ii;
            
        elseif pentaflag==1
            %% PENTA 
            mappasl=repmat(elemconv(ii,1:6),1,3);
            mappasl=mappasl(1,doforder6)*3+repmat([-2:1:0],1,6);
            map(ii,1:18)=mappasl;
            
            map_penta(numpenta,:)=mappasl;
            listpenta(numpenta,1)=ii;
            
        end
        
    end % Main element loop end
    
    if numhexa ~= Hexanum
        error('FATAL')
    end
    
    %% HEXA assembly
    row=reshape(repmat(map_hexa,1,24)',24,24,numhexa);
    col=permute(row,[2 1 3]);
    KassH=sparse(col(:),row(:),Ktemp(:));
    MassH=sparse(col(:),row(:),Mtemp(:));
    MasslH=diag(sum(MassH));
    
    if GyroFlag==1
        GyroH=sparse(col(:),row(:),Gtemp(:));
        CorioH=sparse(col(:),row(:),Ctemp(:));
    end
    
    if numpenta >=1
        if numpenta ~= Pentanum
            error('FATAL')
        end
        
        % only penta mesh
        if Hexanum == 0
            KassH=sparse(nnodi*3,nnodi*3,0);
            MassH=sparse(nnodi*3,nnodi*3,0);
            MasslH=sparse(nnodi*3,nnodi*3,0);
            if GyroFlag==1
                GyroH=sparse(nnodi*3,nnodi*3,0);
                CorioH=sparse(nnodi*3,nnodi*3,0);
            end
        end
        
        NM=max(max(max(row)));
        if NM<=nnodi*3-1
            % expand the matrix to fit the whole dimension if needed
            % can happen only if there are penta elements, not for full hexa
            % mesh
            Kt=sparse(nnodi*3,nnodi*3);
            Mt=sparse(nnodi*3,nnodi*3);
            Mtl=sparse(nnodi*3,nnodi*3);
            if GyroFlag==1
                Gt=sparse(nnodi*3,nnodi*3);
                Ct=sparse(nnodi*3,nnodi*3);
                Gt(1:NM,1:NM)=GyroH;
                Ct(1:NM,1:NM)=CorioH;
                clear GyroH CorioH
                GyroH=Gt;
                CorioH=Ct;
                clear Gt Ct
            end
            Kt(1:NM,1:NM)=KassH;
            Mt(1:NM,1:NM)=MassH;
            Mtl(1:NM,1:NM)=MasslH;
            clear KassH MassH MasslH
            KassH=Kt;
            MassH=Mt;
            MasslH=Mtl;
            clear Kt Mt Mtl row col
        end
        
        %% PENTA assembly
        row=reshape(repmat(map_penta,1,18)',18,18,numpenta);
        col=permute(row,[2 1 3]);
        KassP=sparse(col(:),row(:),KtempP(:));
        MassP=sparse(col(:),row(:),MtempP(:));
        MasslP=diag(sum(MassP));
        
        if GyroFlag==1
            GyroP=sparse(col(:),row(:),GtempP(:));
            CorioP=sparse(col(:),row(:),CtempP(:));
        end
        
        NM=max(max(max(row)));
        if NM<=nnodi*3-1
            % expand the matrix to fit the whole dimension if needed
            % can happen only if there are penta elements, not for full
            % hexa mesh
            Kt=sparse(nnodi*3,nnodi*3);
            Mt=sparse(nnodi*3,nnodi*3);
            Mtl=sparse(nnodi*3,nnodi*3);
            
            if GyroFlag==1
                Gt=sparse(nnodi*3,nnodi*3);
                Ct=sparse(nnodi*3,nnodi*3);
                Gt(1:NM,1:NM)=GyroP;
                Ct(1:NM,1:NM)=CorioP;
                clear GyroP CorioP
                GyroP=Gt;
                CorioP=Ct;
                clear Gt Ct
            end
            
            Kt(1:NM,1:NM)=KassP;
            Mt(1:NM,1:NM)=MassP;
            Mtl(1:NM,1:NM)=MasslP;
            clear KassP MassP MasslP
            KassP=Kt;
            MassP=Mt;
            MasslP=Mtl;
            clear Kt Mt Mtl row col
        end
        
        Kass=KassH+KassP;
        Mass=MassH+MassP;
        Massl=MasslH+MasslP;
        if GyroFlag==1
            Gyro= GyroH+GyroP;
            Corio= CorioH+CorioP;
        end
        
    else
        Kass=KassH;
        Mass=MassH;
        Massl=MasslH;
        
        if GyroFlag==1
            Gyro= GyroH;
            Corio= CorioH;
        end
    end
    
    % fine parte di programma per elementi del primo ordine %%%%%%
    
elseif order ==2
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% SECOND ORDER ELEMENTS ----------------------------------------------------
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % creazione mappa elementi con numerazione nodi matlab
    
    IDnodi=nodi2(:,1);
    elemconv=zeros(nele,20);
    distorsione=zeros(nele,2);
    numhexa=0;
    numpenta=0;
    
    if GyroFlag==0
        Gyro= 'Not Requested';
        Corio='Not Requested';
    end
    
    % preallocate matrices
    % Strain displacement matrix at element centroid
    BK=zeros(6,60,nele);
    % Stress displacement matrices at gauss point locations
    BG=zeros(120,60,nele);
    SG=zeros(120,60,nele);
    map=zeros(nele,60);
    
    % transformation matrices dummy per output
    Te=0;
    
    % Volume
    % vtot = [ Total Volume (hexa+penta), Penta only volume];
    vtot=[0 0];
    
    % temporary matrices
    [Hexanum,Pentanum,ListPenta]=GetElemTypes(elementi2);
    % penta
    KtempP=zeros(45,45,Pentanum);
    MtempP=zeros(45,45,Pentanum);
    
    [Ns21,Nt21,Nu21]= ShapeFunsP15;
    
    if GyroFlag==1
        GtempP=zeros(45,45,Pentanum);
        CtempP=zeros(45,45,Pentanum);
    end
    map_penta=zeros(Pentanum,45);
    listpenta=zeros(Pentanum,1);
    
    %hexa
    Ktemp=zeros(60,60,Hexanum);
    Mtemp=zeros(60,60,Hexanum);
    
    [n27,Ns27,Nt27,Nu27,Ns8,Nt8,Nu8]= ShapeFunsH20;
    
    if GyroFlag==1
        Gtemp=zeros(60,60,Hexanum);
        Ctemp=zeros(60,60,Hexanum);
    end
    map_hexa=zeros(Hexanum,60);
    listhexa=zeros(Hexanum,1);
    
    doforder20=[1:20:60 2:20:60 3:20:60 4:20:60 5:20:60 6:20:60 7:20:60 8:20:60 9:20:60 10:20:60  ...
        11:20:60 12:20:60 13:20:60 14:20:60 15:20:60 16:20:60 17:20:60 18:20:60 19:20:60 20:20:60];
    doforder15=[1:15:45 2:15:45 3:15:45 4:15:45 5:15:45 6:15:45 7:15:45 ...
        8:15:45 9:15:45 10:15:45 11:15:45 12:15:45 13:15:45 14:15:45 15:15:45];
    
    [~,elemconv(:,:)] = builtin('_ismemberhelper',elementi(:,2:21),nodi2(:,1));
    
    % Element assembly loop
    for ii=1:nele

%         [~,elemconv(ii,:)] = ismember(elementi(ii,2:21),nodi2(:,1));

        % genera la matrice delle coordinate locali pescando dalla matrice
        % delle coordinate dei nodi
        coord(1:20,1:3)=nodi2(elemconv(ii,:),2:4);
        
        pentaflag=0; % 0 Hexa 1 Penta
        
        if elemconv(ii,3) == elemconv(ii,4) && elemconv(ii,7) == elemconv(ii,8)
            % penta element
            pentaflag=1; % is penta
            % adjust elemconv table
            elemconv(ii,1:20)=[elemconv(ii,[1:3,5:7,9:10,12:14,16:19]), 0,0,0,0,0];
            
            % adjust coordinates
            coord15=coord([1:3,5:7,9:10,12:14,16:19],:);
            
            if multimaterial == 1
                Ei=E(elem_material(ii,2));
                NUi=NU(elem_material(ii,2));
                rhoi=rho(elem_material(ii,2));
            end
            
            [Kel,Melc,Gel,Cel,Vel,DPel,btotg,B,Delm]=Penta15(Ei,NUi,rhoi,coord15,n27,Ns27,Nt27,Nu27,Ns21,Nt21,Nu21,GyroFlag);
            numpenta=numpenta+1;
            KtempP(:,:,numpenta)=Kel;
            MtempP(:,:,numpenta)=Melc;
            
            if GyroFlag==1
                GtempP(:,:,numpenta)=Gel;
                CtempP(:,:,numpenta)=Cel;
            end
            % strain-displacement matrix at centroid
            BK(1:6,1:45,ii)=B;
            % stress-displacement matrix at integration locations
            SG(1:6,1:45,ii)=Delm*btotg(:,:,1);
            SG(7:12,1:45,ii)=Delm*btotg(:,:,2);
            SG(13:18,1:45,ii)=Delm*btotg(:,:,3);
            SG(19:24,1:45,ii)=Delm*btotg(:,:,4);
            SG(25:30,1:45,ii)=Delm*btotg(:,:,5);
            SG(31:36,1:45,ii)=Delm*btotg(:,:,6);
            SG(37:42,1:45,ii)=Delm*btotg(:,:,7);
            SG(43:48,1:45,ii)=Delm*btotg(:,:,8);
            SG(49:54,1:45,ii)=Delm*btotg(:,:,9);
            SG(55:60,1:45,ii)=Delm*btotg(:,:,10);
            SG(61:66,1:45,ii)=Delm*btotg(:,:,11);
            SG(67:72,1:45,ii)=Delm*btotg(:,:,12);
            SG(73:78,1:45,ii)=Delm*btotg(:,:,13);
            SG(79:84,1:45,ii)=Delm*btotg(:,:,14);
            SG(85:90,1:45,ii)=Delm*btotg(:,:,15);
            
            vtot=vtot+[Vel, Vel];
            distorsione(ii,:)=[elementi(ii,1),DPel];
            
        elseif elemconv(ii,3) ~= elemconv(ii,4) && elemconv(ii,7) ~= elemconv(ii,8) % hexa element
            
            if multimaterial == 1
                Ei=E(elem_material(ii,2));
                NUi=NU(elem_material(ii,2));
                rhoi=rho(elem_material(ii,2));
            end
            
            % calcolo matrici dell'elemento
            [Kel,Melc,Gel,Cel,Vel,DPel,btotg,B,Delm]=Hexa20(Ei,NUi,rhoi,coord,n27,Ns27,Nt27,Nu27,Ns8,Nt8,Nu8,GyroFlag);
            numhexa=numhexa+1;
            Ktemp(:,:,numhexa)=Kel;
            Mtemp(:,:,numhexa)=Melc;

            if GyroFlag==1
                Gtemp(:,:,numhexa)=Gel;
                Ctemp(:,:,numhexa)=Cel;
            end
            
            % strain-displacement matrix at centroid
            BK(:,:,ii)=B;
            % stress-displacement matrix at integration locations
            SG(1:6,:,ii)=Delm*btotg(:,:,1);
            SG(7:12,:,ii)=Delm*btotg(:,:,2);
            SG(13:18,:,ii)=Delm*btotg(:,:,3);
            SG(19:24,:,ii)=Delm*btotg(:,:,4);
            SG(25:30,:,ii)=Delm*btotg(:,:,5);
            SG(31:36,:,ii)=Delm*btotg(:,:,6);
            SG(37:42,:,ii)=Delm*btotg(:,:,7);
            SG(43:48,:,ii)=Delm*btotg(:,:,8);
            %             SG(49:54,:,ii)=Delm*btotg(:,:,9);
            %             SG(55:60,:,ii)=Delm*btotg(:,:,10);
            %             SG(61:66,:,ii)=Delm*btotg(:,:,11);
            %             SG(67:72,:,ii)=Delm*btotg(:,:,12);
            %             SG(73:78,:,ii)=Delm*btotg(:,:,13);
            %             SG(79:84,:,ii)=Delm*btotg(:,:,14);
            %             SG(85:90,:,ii)=Delm*btotg(:,:,15);
            %             SG(91:96,:,ii)=Delm*btotg(:,:,16);
            %             SG(97:102,:,ii)=Delm*btotg(:,:,17);
            %             SG(103:108,:,ii)=Delm*btotg(:,:,18);
            %             SG(109:114,:,ii)=Delm*btotg(:,:,19);
            %             SG(115:120,:,ii)=Delm*btotg(:,:,20);
            %
            vtot(1)=vtot(1)+Vel;
            distorsione(ii,:)=[elementi(ii,1),DPel];
            
        else
            Numerr=num2str(ii);
            error(['Element',' ', Numerr,' contains a number of nodes that does ', ...
                ' ','not correspond to an Hexahedron element (20 nodes) or a Pentahedron', ...
                ' element (15 nodes). Please review your mesh.']);
            
        end
        
        if pentaflag==0
            %% Hexa Assembly
            % 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,doforder20)*3+repmat([-2:1:0],1,20);
            map(ii,:)=mappasl;
            
            map_hexa(numhexa,:)=mappasl;
            listhexa(numhexa,1)=ii;
            
        elseif pentaflag==1
            %% PENTA Assembly
            mappasl=repmat(elemconv(ii,1:15),1,3);
            mappasl=mappasl(1,doforder15)*3+repmat([-2:1:0],1,15);
            map(ii,1:45)=mappasl;
            
            map_penta(numpenta,:)=mappasl;
            listpenta(numpenta,1)=ii;
            
        end
        
    end % Main element loop end
    
    if numhexa ~= Hexanum
        error('FATAL')
    end
    %% HEXA assembly
    row=reshape(repmat(map_hexa,1,60)',60,60,numhexa);
    col=permute(row,[2 1 3]);
    KassH=sparse(col(:),row(:),Ktemp(:));
    MassH=sparse(col(:),row(:),Mtemp(:));
    MasslH=diag(sum(MassH));
    
    if GyroFlag==1
        GyroH=sparse(col(:),row(:),Gtemp(:));
        CorioH=sparse(col(:),row(:),Ctemp(:));
    end
    
    if numpenta >=1
        
        % only penta mesh
        if Hexanum == 0
            KassH=sparse(nnodi*3,nnodi*3,0);
            MassH=sparse(nnodi*3,nnodi*3,0);
            MasslH=sparse(nnodi*3,nnodi*3,0);
            if GyroFlag==1
                GyroH=sparse(nnodi*3,nnodi*3,0);
                CorioH=sparse(nnodi*3,nnodi*3,0);
            end
        end
        
        if numpenta ~= Pentanum
            error('FATAL')
        end
        NM=max(max(max(row)));
        if NM<=nnodi*3-1
            % expand the matrix to fit the whole dimension if needed
            % can happen only if there are penta elements, not for full hexa
            % mesh
            Kt=sparse(nnodi*3,nnodi*3);
            Mt=sparse(nnodi*3,nnodi*3);
            Mtl=sparse(nnodi*3,nnodi*3);
            if GyroFlag==1
                Gt=sparse(nnodi*3,nnodi*3);
                Ct=sparse(nnodi*3,nnodi*3);
                Gt(1:NM,1:NM)=GyroH;
                Ct(1:NM,1:NM)=CorioH;
                clear GyroH CorioH
                GyroH=Gt;
                CorioH=Ct;
                clear Gt Ct
            end
            Kt(1:NM,1:NM)=KassH;
            Mt(1:NM,1:NM)=MassH;
            Mtl(1:NM,1:NM)=MasslH;
            clear KassH MassH MasslH
            KassH=Kt;
            MassH=Mt;
            MasslH=Mtl;
            clear Kt Mt Mtl row col
        end
        
        %% PENTA assembly
        row=reshape(repmat(map_penta,1,45)',45,45,numpenta);
        col=permute(row,[2 1 3]);
        KassP=sparse(col(:),row(:),KtempP(:));
        MassP=sparse(col(:),row(:),MtempP(:));
        MasslP=diag(sum(MassP));
        
        if GyroFlag==1
            GyroP=sparse(col(:),row(:),GtempP(:));
            CorioP=sparse(col(:),row(:),CtempP(:));
        end
        
        NM=max(max(max(row)));
        if NM<=nnodi*3-1
            % expand the matrix to fit the whole dimension if needed
            % can happen only if there are penta elements, not for full
            % hexa mesh
            Kt=sparse(nnodi*3,nnodi*3);
            Mt=sparse(nnodi*3,nnodi*3);
            Mtl=sparse(nnodi*3,nnodi*3);
            
            if GyroFlag==1
                Gt=sparse(nnodi*3,nnodi*3);
                Ct=sparse(nnodi*3,nnodi*3);
                Gt(1:NM,1:NM)=GyroP;
                Ct(1:NM,1:NM)=CorioP;
                clear GyroP CorioP
                GyroP=Gt;
                CorioP=Ct;
                clear Gt Ct
            end
            
            Kt(1:NM,1:NM)=KassP;
            Mt(1:NM,1:NM)=MassP;
            Mtl(1:NM,1:NM)=MasslP;
            clear KassP MassP MasslP
            KassP=Kt;
            MassP=Mt;
            MasslP=Mtl;
            clear Kt Mt Mtl row col
        end
        
        Kass=KassH+KassP;
        Mass=MassH+MassP;
        Massl=MasslH+MasslP;
        if GyroFlag==1
            Gyro= GyroH+GyroP;
            Corio= CorioH+CorioP;
        end
        
    else
        Kass=KassH;
        Mass=MassH;
        Massl=MasslH;
        
        if GyroFlag==1
            Gyro= GyroH;
            Corio= CorioH;
        end
    end
end % fine parte di programma per elementi del secondo ordine %%%%%%%%%%%%%%%


end % poliFEMo end
