%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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[p1_neighbours, p2_neighbours, p1b_neighbours] = find_struct_leftneighbour(MESH, PARAMETERS)

	cprec = 1e-8;

	if (max(strcmp(PARAMETERS.DOMAIN.GEOMETRY, {'RECTANGLE', 'DIHEDRON', 'STEP'})) ...
					&& strcmp(PARAMETERS.MESH.GENERATION, 'NS2DDV'))
		switch PARAMETERS.DOMAIN.GEOMETRY
			case {'DIHEDRON'}
				dx = 1./PARAMETERS.MESH.NBSEG_PERUNIT_X;
				npy_p1 = PARAMETERS.MESH.NBNODES_INBL_Y + PARAMETERS.MESH.NBSEG_OUTBL_Y;
				npy_p2 = 2*npy_p1-1;
				npy_p1b = npy_p1;
			case {'RECTANGLE'}
				dx = (PARAMETERS.DOMAIN.XMAX-PARAMETERS.DOMAIN.XMIN)/PARAMETERS.MESH.NBSEG_X;
				npy_p1 = PARAMETERS.MESH.NBSEG_Y+1;
				npy_p2 = 2*npy_p1-1;
				npy_p1b = npy_p1;
			case {'STEP'}
				error('NS2DDV mesh generation is not implemented yet for STEP geometry');
			otherwise
				error('Problem in find_struct_leftneighbour');
		end
	else
		error('Problem in find_struct_leftneighbour');
	end
	Lx = PARAMETERS.DOMAIN.XMAX-PARAMETERS.DOMAIN.XMIN;

	nbneighbours = 6;

	% P2-nodes on {xmax} x [ymin,ymax]
	p = MESH.p2;
	np = MESH.np2;
	npy = npy_p2;
	indpref = zeros(npy,1);
	k = 1;
	for i=1:np
		ploc = p(:,i);
		if (abs(ploc(1)-PARAMETERS.DOMAIN.XMAX)/Lx < cprec)
			indpref(k) = i;
			k = k+1;
		end
	end

	if ~(k-1 == npy)
		error('Problem in find_struct_leftneighbour: it seems that the geometry or the mesh is not well suited for running this routine')
	end
	p_neighbours = cell(nbneighbours+1,1);
	p_neighbours{1} = indpref;
	for kn=1:nbneighbours
		p_neighbours{kn+1} = zeros(np,1);
		for k=1:npy
			ploc_neighbour_candidate = p(:,indpref(k)) - [kn*dx; 0.];
			indpneighbour = 0;
			for i=1:np
				if (sqrt((ploc_neighbour_candidate(1)-p(1,i))^2+(ploc_neighbour_candidate(2)-p(2,i))^2) < cprec)
					indpneighbour = i;
					break;
				end
			end
			if indpneighbour
				p_neighbours{kn+1}(indpref(k),1) = indpneighbour;
			else
				error('Problem in find_struct_leftneighbour: it seems that the geometry or the mesh is not well suited for running this routine');
			end
		end
	end
	p2_neighbours = p_neighbours;

	% P1-nodes on {xmax} x [ymin,ymax]
	p = MESH.p1;
	np = MESH.np1;
	npy = npy_p1;
	indpref = zeros(npy,1);
	k = 1;
	for i=1:np
		ploc = p(:,i);
		if (abs(ploc(1)-PARAMETERS.DOMAIN.XMAX)/Lx < cprec)
			indpref(k) = i;
			k = k+1;
		end
	end
	if ~(k-1 == npy)
		error('Problem in find_struct_leftneighbour: it seems that the geometry or the mesh is not well suited for running this routine');
	end
	p_neighbours = cell(nbneighbours+1,1);
	p_neighbours{1} = indpref;
	for kn=1:nbneighbours
		p_neighbours{kn+1} = zeros(np,1);
		for k=1:npy
			ploc_neighbour_candidate = p(:,indpref(k)) - [kn*dx; 0.];
			indpneighbour = 0;
			for i=1:np
				if (sqrt((ploc_neighbour_candidate(1)-p(1,i))^2+(ploc_neighbour_candidate(2)-p(2,i))^2) < cprec)
					indpneighbour = i;
					break;
				end
			end
			if indpneighbour
				p_neighbours{kn+1}(indpref(k),1) = indpneighbour;
			else
				error('Problem in find_struct_leftneighbour: it seems that the geometry or the mesh is not well suited for running this routine');
			end
		end
	end
	p1_neighbours = p_neighbours;

	% P1b-nodes on {xmax} x [ymin,ymax]
	p = MESH.p1b;
	np = MESH.np1b;
	npy = npy_p1b;
	indpref = zeros(npy,1);
	k = 1;
	for i=1:np
		ploc = p(:,i);
		if (abs(ploc(1)-PARAMETERS.DOMAIN.XMAX)/Lx < cprec)
			indpref(k) = i;
			indpref(k) = i;
			k = k+1;
		end
	end
	if ~(k-1 == npy)
		error('Problem in find_struct_leftneighbour: it seems that the geometry or the mesh is not well suited for running this routine');
	end
	p_neighbours = cell(nbneighbours+1,1);
	p_neighbours{1} = indpref;
	for kn=1:nbneighbours
		p_neighbours{kn+1} = zeros(np,1);
		for k=1:npy
			ploc_neighbour_candidate = p(:,indpref(k)) - [kn*dx; 0.];
			indpneighbour = 0;
			for i=1:np
				if (sqrt((ploc_neighbour_candidate(1)-p(1,i))^2+(ploc_neighbour_candidate(2)-p(2,i))^2) < cprec)
					indpneighbour = i;
					break;
				end
			end
			if indpneighbour
				p_neighbours{kn+1}(indpref(k),1) = indpneighbour;
			else
				error('Problem in find_struct_leftneighbour: it seems that the geometry or the mesh is not well suited for running this routine');
			end
		end
	end
	p1b_neighbours = p_neighbours;

end