%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of NS2DDV.                                                      %
%                                                                                   %
% Copyright(C) 2011-2018    C. Calgaro  (caterina.calgaro@math.univ-lille1.fr)      %
%                           E. Creusé   (emmanuel.creuse@math.univ-lille1.fr)       %
%                           T. Goudon   (thierry.goudon@inria.fr)                   %
%                           A. Mouton   (alexandre.mouton@math.univ-lille1.fr)      %
%                                                                                   %
% NS2DDV is free software: you can redistribute it and/or modify it under the terms %
% of the GNU General Public License as published by the Free Software Foundation,   %
% either version 3 of the License, or (at your option) any later version.           %
%                                                                                   %
% NS2DDV is distributed in the hope that it will be useful, but WITHOUT ANY         %
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   %
% PARTICULAR PURPOSE. See the GNU General Public License for more details.          %
%                                                                                   %
% You should have received a copy of the GNU General Public License along with      %
% NS2DDV. If not, see <http://www.gnu.org/licenses/>.                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[p,t] = structmesh_step_2d(PARAMETERS)

	% INPUT
	% PARAMETERS	List of simulation parameters
	%
	% OUTPUT
	% p 			P1 nodes
	% t				Triangles described by P1 nodes

	dx = (PARAMETERS.DOMAIN.XMAX-PARAMETERS.DOMAIN.XMIN)/PARAMETERS.MESH.NBSEG_X;
	dy = (PARAMETERS.DOMAIN.YMAX-PARAMETERS.DOMAIN.YMIN)/PARAMETERS.MESH.NBSEG_Y;

	Nx = floor((PARAMETERS.DOMAIN.XSTEP-PARAMETERS.DOMAIN.XMIN)/dx);
	Ny = floor((PARAMETERS.DOMAIN.YSTEP-PARAMETERS.DOMAIN.YMIN)/dy);

	ddx = (PARAMETERS.DOMAIN.XSTEP-PARAMETERS.DOMAIN.XMIN)/Nx;
	ddy = (PARAMETERS.DOMAIN.YSTEP-PARAMETERS.DOMAIN.YMIN)/Ny;

	dddx = (PARAMETERS.DOMAIN.XMAX-PARAMETERS.DOMAIN.XSTEP)/Nx;
	dddy = (PARAMETERS.DOMAIN.YMAX-PARAMETERS.DOMAIN.YSTEP)/Ny;

	ind = -ones(PARAMETERS.MESH.NBSEG_X+1,PARAMETERS.MESH.NBSEG_Y+1);

	np = (PARAMETERS.MESH.NBSEG_X+1)*(PARAMETERS.MESH.NBSEG_Y-Ny+1) + ...
					(PARAMETERS.MESH.NBSEG_X-Nx+1)*PARAMETERS.MESH.NBSEG_Y;
	p = zeros(2,np);

	switch PARAMETERS.DOMAIN.STEP_POSITION
		case {'LEFT'}
			% Upper left part of the domain
			k = 1;
			for i=0:Nx
				for j=Ny:PARAMETERS.MESH.NBSEG_Y
					ind(i+1,j+1) = k;
					p(1,k) = PARAMETERS.DOMAIN.XMIN + i*ddx;
					p(2,k) = PARAMETERS.DOMAIN.YSTEP + (j-Ny)*dddy;
					k = k+1;
				end
			end
			% Lower right part of the domain
			for i=Nx:PARAMETERS.MESH.NBSEG_X
				for j=0:(Ny-1)
					ind(i+1,j+1) = k;
					p(1,k) = PARAMETERS.DOMAIN.XSTEP + (i-Nx)*dddx;
					p(2,k) = PARAMETERS.DOMAIN.YMIN + j*ddy;
					k = k+1;
				end
			end
			% Upper right part of the domain
			for i=(Nx+1):PARAMETERS.MESH.NBSEG_X
				for j=Ny:PARAMETERS.MESH.NBSEG_Y
					ind(i+1,j+1) = k;
					p(1,k) = PARAMETERS.DOMAIN.XSTEP + (i-Nx)*dddx;
					p(2,k) = PARAMETERS.DOMAIN.YSTEP + (j-Ny)*dddy;
					k = k+1;
				end
			end

			disp(k-1)
			disp(np)

		case {'RIGHT'}



		otherwise
			error('Unknown position of step in STEP domain geometry');
	end	

end