function msh=mesh_generator(geo,siz,d,msh)
nn0=msh.nn;  % eddigi connectivity matrix
xy0=msh.xy;  % eddigi koordinatamatrix

nop=size(geo,2); % partok szama

%% ÚJ PARTOK KOORDINATA ÉS CONNECTIVITY MATRIXANAK ELOALLITASA %%%%%%%%%%%%%%%%%

for p=1:nop % for ciklus a partokra
  
  % 1. aktualis part koordinatamatrixanak eloallitasa:
  
  % sarokpontok koordinatai
  A=geo{p}(1,:);      % A----------B
  B=geo{p}(2,:);      % |          |
  C=geo{p}(3,:);      % |          |
  D=geo{p}(4,:);      % C----------D
  r=siz{p}(1);        % sorok szama
  c=siz{p}(2);        % oszlopok szama
  
  noe=r*d;                 % elemek szama
  non1=d*c+1;              % csomopontok szama vizszintes vonalakon
  non2=(d-1)*(c+1);        % koztes csomopontok szama
  non=(r+1)*non1+r*non2;   % csomopontok szama
  
  idx=[1:d*r+1]';          % összes fuggoleges pont index
  idx1=[1:d:d*r+1]';       % sarokpont index
  idx2=setdiff(idx,idx1);  % belsopont index
  
  ACxy=A+(C-A).*[0:r*d]'./(r*d); % koordinatak bal oldalon
  BDxy=B+(D-B).*[0:r*d]'./(r*d); % koordinatak jobb oldalon
  ACxy1=ACxy(idx1,:);            % sarok koordinatak bal oldalon
  BDxy1=BDxy(idx1,:);            % sarok koordinatak jobboldalon
  ACxy2=ACxy(idx2,:);            % belso koordinatak baloldalon
  BDxy2=BDxy(idx2,:);            % belso koordinatak jobboldalon
  
  ACxy1=repelem(ACxy1,non1,1);   % bal sarok koordinatak ismetlese non1-szer
  BDxy1=repelem(BDxy1,non1,1);   % jobb sarok koordinatak ismetlese non1-szer
  ACxy2=repelem(ACxy2,c+1,1);    % bal belso koordinatak ismetlese oszlop+1-szer
  BDxy2=repelem(BDxy2,c+1,1);    %jobb belso koordinatak ismetlese oszlop+1-szer
  
  idx=[1:non]';                                     % index 1:csomopontszamig
  idx1=[1:(non1+non2):r*(non1+non2)+1]+[0:c*d]';    % indexek vizszintes vonalakon
  idx2=setdiff(idx,idx1);                           % belso indexek
  
  % az aktualis koordinamatrix:
  xy1{p,1}(idx1,:)=ACxy1+(BDxy1-ACxy1).*repmat([0:c*d]'./(c*d),r+1,1);
  xy1{p,1}(idx2,:)=ACxy2+(BDxy2-ACxy2).*repmat([0:c]'./c,(d-1)*r,1);
  
  % 2. az aktualis part connectivity matrixanak eloallitasa:
  for i=1:r
    for j=1:c
      for k=1:d+1
        nn1{p,1}((i-1)*c+j,k)=i*(non1+non2)+(j-1)*d+k;
        nn1{p,1}((i-1)*c+j,3*d+2-k)=(i-1)*(non1+non2)+(j-1)*d+k;
      end
      for k=1:d-1
        nn1{p,1}((i-1)*c+j,2*d+1-k)=i*non1+(i-1)*non2+(k-1)*(c+1)+j+1;
        nn1{p,1}((i-1)*c+j,3*d+1+k)=i*non1+(i-1)*non2+(k-1)*(c+1)+j;
      end
    end
   end
   
end % for p

%% PARTOK OSSZERAGASZTASA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;

% 1. kozos koordinatamatrix letrehozasa:
xyA=[xy0;cell2mat(xy1)];         % osszes koordinamatrix egymas ala rakva
if isOctave==0
    [xy,ia,ic]=uniquetol(xyA,'ByRows',true);   % egyezo koordinatak torlese
else
    xyA=round(xyA*1000)/1000;
    [xy,ia,ic]=unique(xyA,'rows');   % egyezo koordinatak torlese
end

ia2=sort(ia);                    % xyzA egyedi sorainak a sorszama
xy=xyA(ia2,:);                   % a kozos koordinatamatrix

% 2. kozos connectivity matrix letrehozasa:
if isOctave==0
    [~,ic2]=ismembertol(xyA,xy,'ByRows',true); % xyzA sorai, hova keruljenek xyz-ben
else
    [~,ic2]=ismember(xyA,xy,'rows'); % xyzA sorai, hova keruljenek xyz-ben
end
sn=size(xy0,1)+1;                % part elso pontjanak inicializalasa
for i=1:nop                      % for ciklus az uj partokra
    nonpp=size(xy1{i},1);        % aktualis part csomopontjainak szama
    newnum=ic2(sn:sn+nonpp-1);   % az aktualis part uj sorszamai
    nnn{i,1}=zeros(size(nn1{i}));% aktualis part nn matrixanak inicializalasa
    for j=1:nonpp                % for ciklus az aktualis part osszes pontjara
        nnn{i,1}(nn1{i}==j)=newnum(j);% aktualis part feltoltese uj sorszamokkal
    end
    sn=sn+nonpp;                 % kovetkezo part kezdosorszama
end

% matrixok msh strukturaba rakasa
msh.xy=xy;
msh.nn=[nn0;cell2mat(nnn)];

% halo abrazolasa
figure
set(gcf,'color','w');
patch('faces',msh.nn,'vertices',msh.xy,'facecolor',[0.7,0.7,0.7],'edgecolor','k')
axis equal
xlim([min(xy(:,1)) max(xy(:,1))]);
ylim([min(xy(:,2)) max(xy(:,2))]);