%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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[boundary_edges,boundary_edges_cells,boundary_edges_length,boundary_edges_normals] = mesh_stars_fv_part2(MESH)

	% INPUT
	% MESH						An incomplete mesh metadata
	%
	% OUTPUT
	% boundary_edges			Array which gives global edges number on the boundary
	% boundary_edges_cells		Contains the number of the triangle which contains a border edge
	% boundary_edges_length		Length of bound edges
	% boundary_edges_normals	Normal unit vector to the bound edges

	nb_boun_edges = 0;
	for ia = 1:size(MESH.ne1full,2)
		if (MESH.e1full(3,ia) ~= 0)
			nb_boun_edges = nb_boun_edges+1;
		end
	end

	boundary_edges        = zeros(1,nb_boun_edges);
	boundary_edges_cells  = zeros(1,nb_boun_edges);
	boundary_edges_length = zeros(1,nb_boun_edges);
	boundary_edges_normals= zeros(2,nb_boun_edges);

	nb_ar_bord=0;
	for ke = 1:MESH.nt1
		for j = 1:3
			ia = abs(MESH.t1e(j,ke));
			if (MESH.e1full(3,ia) ~= 0)
				nb_ar_bord = nb_ar_bord+1;
				boundary_edges(1,nb_ar_bord) = ia;
				boundary_edges_cells(1,nb_ar_bord) = ke;

				jp1 = mod(j  ,3)+1;
				jm1 = mod(jp1,3)+1;
				neup = MESH.t1(jp1,ke);
				neum = MESH.t1(jm1,ke);

				normal_bord_1 = -(MESH.p1(2,neup) - MESH.p1(2,neum));
				normal_bord_2 =  (MESH.p1(1,neup) - MESH.p1(1,neum));
				boundary_edges_length(1,nb_ar_bord) = sqrt(normal_bord_1^2+normal_bord_2^2);

				vect_int_1 = (MESH.p1(1,MESH.t1(j,ke)) - MESH.p1(1,neup));
				vect_int_2 = (MESH.p1(2,MESH.t1(j,ke)) - MESH.p1(2,neup));
            
				if normal_bord_1*vect_int_1+normal_bord_2*vect_int_2 > 0
					boundary_edges_normals(1,nb_ar_bord) = -normal_bord_1;
					boundary_edges_normals(2,nb_ar_bord) = -normal_bord_2;
				else
					boundary_edges_normals(1,nb_ar_bord) =  normal_bord_1;
					boundary_edges_normals(2,nb_ar_bord) =  normal_bord_2;
				end
			end
		end
	end

end