function [nodi2,elementi2]=RemeshTo2ndOrderElements(nodi,elementi)

% vettore connessioni hexa - 12 edge
% [nodo 1, nodo 2 , nodo addizionale]
hexav=[1 2 9;
    2 3 10;
    3 4 11;
    4 1 12;
    5 6 13;
    6 7 14;
    7 8 15;
    8 5 16;
    1 5 17;
    2 6 18;
    3 7 19;
    4 8 20];

% vettori connessioni penta - 9 edge
% [nodo 1, nodo 2 , nodo addizionale]
pentav=[ 1 2 7;
    2 3 8;
    3 1 9;
    4 5 10;
    5 6 11;
    6 4 12;
    1 4 13;
    2 5 14;
    3 6 15];

nele=size(elementi,1);
elementi2=zeros(nele,21);
elementi2(:,1:9)=elementi;
nnodi=size(nodi,1);
% numero di nodi non modificati
nnodi_or=nnodi;
% numero del nodo più alto
num_nod=max(nodi(:,1));
% contatore di inizio per numerare i nodi aggiunti
num_cont=num_nod+1;
contat=1; % contatore per loop sugli edge
% vettore dei nodi aggiunti
nodi_add=[0 0 0 0];
empty=nodi_add;

% inizio loop di split degli edge di tutti gli elementi
for ce=1:nele
    
    elp=elementi(ce,2:9);
    
    if size(unique(elp),2) ==8
        % caso hexa
        numedge=12;
        
        % loop per ogni edge
        for nd=1:numedge
            % calcolo coordinate nodo medio
            xmid=(nodi(elp(1,hexav(nd,1)),2)+nodi(elp(1,hexav(nd,2)),2))/2;
            ymid=(nodi(elp(1,hexav(nd,1)),3)+nodi(elp(1,hexav(nd,2)),3))/2;
            zmid=(nodi(elp(1,hexav(nd,1)),4)+nodi(elp(1,hexav(nd,2)),4))/2;
            % aggiungo le coordinate del punto medio dell'edge
            nodi_add(contat,1:4)= [num_cont, xmid, ymid, zmid];
            % aggiungo il nodo all'elemento
            elementi2(ce,hexav(nd,3))=num_cont;
            
            % aggiorno contatori
            contat=contat+1;
            num_cont=num_cont+1;
        end % fine loop per gli edge
        
        
    elseif size(unique(elp),2) ==7
        % caso penta
        numedge=9;
        
        % loop per ogni edge
        for nd=1:numedge
            % calcolo coordinate nodo medio
            xmid=(nodi(elp(1,pentav(nd,1)),2)+nodi(elp(1,pentav(nd,2)),2))/2;
            ymid=(nodi(elp(1,pentav(nd,1)),3)+nodi(elp(1,pentav(nd,2)),3))/2;
            zmid=(nodi(elp(1,pentav(nd,1)),4)+nodi(elp(1,pentav(nd,2)),4))/2;
            % aggiungo le coordinate del punto medio dell'edge
            nodi_add(contat,1:4)= [num_cont, xmid, ymid, zmid];
            % aggiungo il nodo all'elemento
            elementi2(ce,hexav(nd,3))=num_cont;
            
            % aggiorno contatori
            contat=contat+1;
            num_cont=num_cont+1;
        end % fine loop per gli edge
        
    end

end % fine loop su tutti gli elementi


% tutti gli edge adesso hanno il loro nodo addizionale
% devo rimuovere i duplicati che si creano per gli edge coincidenti tra i
% vari elementi
[nodi_add_tmp,ia,ic] = unique(nodi_add(:,2:4),'rows','stable');

% contatore di inizio per numerare i nodi aggiunti
new_count=num_nod+1;
% nuova numerazione dei nodi aggiunti
nodi_add_new=[(new_count:1:new_count+size(nodi_add_tmp,1)-1)' nodi_add_tmp];



for cl=1:size(ic,1)
   if ic(cl) ~= 0  
    [r1,~]=find(ic(cl) == ic(:));
    if size(r1,1) >=2
        for cl2=1:size(r1,1)
            [r2,c2]=find(nodi_add(r1(cl2),1) == elementi2(:,2:21));
            elementi2(r2,c2)= new_count;
            new_count=new_count+1;
        end
    end
    [r,~]=find(ic(cl) == ic(:));
    ic(r)=0;
   end
end

% % loop per tutti i nodi aggiunti
% for cn=1:size(nodi_add,1)
%     % nodo in esame
%     nod=nodi_add(cn,:);
%     % solo se il nodo non è già stato eliminato da un loop precedente
%     if ~all(nod == 0) %
%         % cerco tutti i nodi che hanno le stesse coordinate del nodo in
%         % esame nod 
%         % nnz(nod(ii,2:4)-nodi_add(:,2:4))~=0
% %         [riga,~]=find(nod(:,2:4) == nodi_add(:,2:4));
%         [riga,~]=find(all(nod(:,2:4)-nodi_add(:,2:4)));
%         % se ci sono corrispondenze
%         if ~isempty(riga)
%         for cn2=1:size(riga,1) % per tutti i nodi trovati corrispondenti
%             % check se il nodo corrispondente trovato non è il nodo stesso
%             % in esame nod
%             if ~(nod(1,1) == nodi_add(riga(cn2),1))
%                 % se non è il nodo in esame 
%                 
%                 %  sostituisco il numero nel nodo dove compare negli
%                 % elementi
%                 [re,ce]=find( nodi_add(riga(cn2),1) == elementi2(:,2:21));
%                 elementi2(re,ce+1) = nod(1,1); % *ones(size(re,1),1) ??
%                 
%                 % e poi lo elimino
%                 nodi_add(riga(cn2),:)=empty;
%                
%             end
%         end % fine loop per nodi corrispondenti
%         end
%     end % fine if
% end % fine loop per tutti i nodi aggiunti   
% 
% % cancello tutti i nodi eliminati e riordino in modo crescente
% nodi_add_tmp=unique(nodi_add);
% nodi_add_tmp(1,:)=[];
% [~,ind]=sort(nodi_add_tmp(:,1),1);
% 
% % aggrego i nodi originali e quelli aggiunti per formare il vettore dei
% % nodi definitivo
% nodi2=[nodi;
%         nodi_add_tmp(ind,:)];
%         
nodi2=[nodi; nodi_add_new];

end %function end