function [xinit,e,ep,ey,ec,ID,IN,Ishape,Ip,Iy,Ic,IDD,IDD1,IDD2] = discret(xl,xt,xpl,xpt,xyl,xyr,xyt,xcl,xcr,xct,h)

% [xinit,e,ep,ey,ec,ID,IN,Ishape,Ip,Iy,Ic,IDD,IDD1,IDD2] = discret(xl,xt,xpl,xpt,xyl,xyr,xyt,xcl,xcr,xct,h)
%          It discretizes a quarter of a domain (due to symmetry)
%          that contains 2 ferromagnetic yokes completed by coils
%          and a domain \Omega_p.
%          x is a vector of nodal coordinates,
%          e contains triples of nodal indices representing triangular
%          elements,
%          ep are indices of the elements within \Omega_p,
%          ey are indices of the elements in the ferromagnetic yoke,
%          ec are indices of the elements within the coil,
%          ID is the Dirichlet boundary, i.e. indices
%          of nodes along \Gamma_D, which is the top, left, and bottom,
%          IN is the Neumann boundary, which is the right part,
%          Idesign contains nodal indices of the design interface,
%          Ip, Iy, and Ic contain nodal indices of the domain
%          where the required magnetic field should be achieved,
%          the domain occupied by the yoke and the one occupied
%          by the coil, respectively.


%sprintf('%s','discret')

% number of design variables
nDesign = 4;

dy = xyt/nDesign;

x1 = xl:h:0;
x2 = 0:h:xt;

epsil = h/10;

nx1 = max(size(x1));
nx2 = max(size(x2));

xinit = [];

e = [];
ep = [];
ey = [];
ec = [];

ID = [];
IN = [];
Ishape = [];
Ip = [];
Iy = [];
Ic = [];
IDD = [];
IDD1 = [];
IDD2 = [];

index = 1;
elindex = 1;

% seznam uzlu a elementu
for i=1:nx1
    for j=1:nx2
        xinit = [xinit [x1(i);x2(j)]];
      
        if i<nx1 & j<nx2
            e = [ e ; [index index+nx2+1 index+1] ; [index index+nx2 index+nx2+1] ];
            elindex = elindex+2;
        end

        if x1(i)>=xpl-epsil & x1(i)<0 & x2(j)<xpt
            ep = [ ep ; elindex-2 ; elindex-1 ];
        end
        
        if x1(i)>=xyl-epsil & x1(i)<xyr & x2(j)<xyt
            ey = [ ey ; elindex-2 ; elindex-1 ];
        end
       
        if x1(i)>=xcl-epsil & x1(i)<xcr & x2(j)>=xyt-epsil & x2(j)<xct
            ec = [ ec ; elindex-2 ; elindex-1 ];
        end

        if j==1 | j==nx2 | i==1
            ID = [ID;index];

            if ~(x1(i)>=xyl-epsil & x1(i)<=xyr+epsil & j==1)
                IDD = [IDD;index];
            else
                IDD2 = [IDD2;index];
            end
        end
        
        if i==nx1
            IN = [IN;index];
            IDD = [IDD;index];
        end
        
        if x1(i)<=xyr+epsil & x1(i)>=xyr-epsil & x2(j)<=xyt+epsil
            Ishape = [Ishape;index];
            IDD = [IDD;index];
        end
        
        if (x1(i)<=xpl+epsil & x1(i)>=xpl-epsil & x2(j)<=xpt+epsil) | (x1(i)>=xpl-epsil & x2(j)<=xpt+epsil & x2(j)>=xpt-epsil)
            Ip = [Ip;index];
            IDD = [IDD;index];
        end
        
        if (x1(i)<=xyl+epsil & x1(i)>=xyl-epsil & x2(j)<=xyt+epsil) | (x1(i)>=xyl-epsil & x1(i)<=xyr+epsil & x2(j)>=xyt-epsil & x2(j)<=xyt+epsil)
            Iy = [Iy;index];
            IDD = [IDD;index];
        end
        
        if (x1(i)<=xcl+epsil & x1(i)>=xcl-epsil & x2(j)>=xyt-epsil & x2(j)<=xct+epsil) | (x1(i)>=xcl-epsil & x1(i)<=xcr+epsil & x2(j)<=xct+epsil & x2(j)>=xct-epsil) | (x1(i)<=xcr+epsil & x1(i)>=xcr-epsil & x2(j)>=xyt-epsil & x2(j)<=xct+epsil)
            Ic = [Ic;index];
            IDD = [IDD;index];
        end
               
        index = index + 1;
    end
end

% pocet uzlu
n = size(xinit,2);

% pocet elementu
ne = size(e,1);
