%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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[e, be, te] = edge_descriptor(t, PARAMETERS)

	% INPUT
	% t 			List of P1 triangles described by P1 nodes
	% PARAMETERS	List of simulation parameters
	%
	% OUTPUT
	% e				Edge descriptor
	% be			Bound edge descriptor
	% te			Triangle descriptor by edges

	% For each edge
	%   e(1:2,i) is the list of P1 nodes that are the extremities
	%   e(3:4,i) are the left and right cells that are separated by the edge
	%   e(5,i)   is the subdomain index in which the edge is included
	%   e(4,i) < 1 if e(:,i) describes a bound edge

	% For each bound edge
	%   be(1:2,i) is the list of P1 nodes that are the extremities
	%   be(3,i)   is the cell that is separated from the exterior by the bound edge
	%   be(4,i)   is dummy
	%   be(5,i)   is the index of the bound edge in the list e above

	% For each triangle
	%  te(k,i)    is the index of the edge in the list e above that is opposed to the P1 node t(k,i) in the triangle

	nt = size(t,2);
	ne = 3*nt;
	e = [];
	be = [];
	te = [];

	% We extract the list of edges from the description of the cells by nodes
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% We consider two temporary edges lists: temp1 will contains the edges with vertex indices (i1,i2) satisfying i1 < i2 and
	% temp2 will contains the edges with vertex indices (i1,i2) satisfying i1 > i2
	% Each edge can appear only one time in temp1
	% Each edge can appear only one time in temp2
	% An interior edge appears one time in temp1 and one time in temp2
	% A bound edge appears one time in the reunion of temp1 and temp2
	% In temp2, we permute i1 and i2 before storing the edge in order to have temp2(*,1) < temp2(*,2)
	% ==> temp1 and temp2 are not completely fullfilled after having read all the cells
	temp1_1 = zeros(nt,5);
	temp1_2 = zeros(nt,5);
	temp1_3 = zeros(nt,5);
	temp2_1 = zeros(nt,5);
	temp2_2 = zeros(nt,5);
	temp2_2 = zeros(nt,5);
	switch PARAMETERS.PARALLELIZATION.TOOLBOX
		case {'SERIAL'}
			for ke=1:nt
				tloc = t(1:3,ke);
				if (tloc(1) < tloc(2))
					temp1_1(ke,:) = [ ...
						tloc(1), ...	% First extremity of the edge
						tloc(2), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(3), ...	% Keep in mind the third vertex of the triangle
						tloc(3)];
				else
					temp2_1(ke,:) = [ ...
						tloc(2), ...	% First extremity of the edge
						tloc(1), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(3), ...	% Keep in mind the third vertex of the triangle
						tloc(3)];
				end
				if (tloc(2) < tloc(3))
					temp1_2(ke,:) = [ ...
						tloc(2), ...	% First extremity of the edge
						tloc(3), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(1), ...	% Keep in mind the third vertex of the triangle
						tloc(1)];
				else
					temp2_2(ke,:) = [ ...
						tloc(3), ...	% First extremity of the edge
						tloc(2), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(1), ...	% Keep in mind the third vertex of the triangle
						tloc(1)];
				end
				if (tloc(3) < tloc(1))
					temp1_3(ke,:) = [ ...
						tloc(3), ...	% First extremity of the edge
						tloc(1), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(2), ...	% Keep in mind the third vertex of the triangle
						tloc(2)];
				else
					temp2_3(ke,:) = [ ...
						tloc(1), ...	% First extremity of the edge
						tloc(3), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(2), ...	% Keep in mind the third vertex of the triangle
						tloc(2)];
				end
			end
		case {'PCT'}
			parfor ke=1:nt
				tloc = t(1:3,ke);
				if (tloc(1) < tloc(2))
					temp1_1(ke,:) = [ ...
						tloc(1), ...	% First extremity of the edge
						tloc(2), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(3), ...	% Keep in mind the third vertex of the triangle
						tloc(3)];
				else
					temp2_1(ke,:) = [ ...
						tloc(2), ...	% First extremity of the edge
						tloc(1), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(3), ...	% Keep in mind the third vertex of the triangle
						tloc(3)];
				end
				if (tloc(2) < tloc(3))
					temp1_2(ke,:) = [ ...
						tloc(2), ...	% First extremity of the edge
						tloc(3), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(1), ...	% Keep in mind the third vertex of the triangle
						tloc(1)];
				else
					temp2_2(ke,:) = [ ...
						tloc(3), ...	% First extremity of the edge
						tloc(2), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(1), ...	% Keep in mind the third vertex of the triangle
						tloc(1)];
				end
				if (tloc(3) < tloc(1))
					temp1_3(ke,:) = [ ...
						tloc(3), ...	% First extremity of the edge
						tloc(1), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(2), ...	% Keep in mind the third vertex of the triangle
						tloc(2)];
				else
					temp2_3(ke,:) = [ ...
						tloc(1), ...	% First extremity of the edge
						tloc(3), ...	% Second extremity of the edge
						ke, ...		% Index of involved triangular cell
						tloc(2), ...	% Keep in mind the third vertex of the triangle
						tloc(2)];
				end
			end
		otherwise
			error(sprintf('Unknown value of PARAMETERS.PARALLELIZATION.TOOLBOX\n'));
	end

	temp1 = [temp1_1; temp1_2; temp1_3];
	newtemp11 = sortrows(temp1);
	for i=1:size(newtemp11,1)
		if (newtemp11(i,1) > 0)
			break;
		end
	end
	newtemp1 = newtemp11(i:size(newtemp11,1),:);

	temp2 = [temp2_1; temp2_2; temp2_3];
	newtemp22 = sortrows(temp2);
	for i=1:size(newtemp22,1)
		if (newtemp22(i,1) > 0)
			break;
		end
	end
	newtemp2 = newtemp22(i:size(newtemp22,1),:);

	clear temp1;
	clear temp1_1;
	clear temp1_2;
	clear temp1_3;
	clear newtemp11;
	clear temp2;
	clear temp2_1;
	clear temp2_2;
	clear temp2_3;
	clear newtemp22;
	

	tempedges = zeros(ne,8);
	n1 = 1;
	n2 = 1;
	%==> The n1 first lines of newtemp1 and the n2 first lines of newtemp2 are made of 0.
	% We skip these lines.
	% We fill the 8 first columns of edges
	i = 1;
	while ((n1 <= size(newtemp1,1)) || (n2 <= size(newtemp2,1)))
		if (n1 > size(newtemp1,1))
			% newtemp1 has been completely scanned
			% The data entirely comes from newtemp2
			tempedges(i,1:3) = newtemp2(n2,1:3);
			tempedges(i,5:6) = newtemp2(n2,4:5);
			n2 = n2+1;
			i = i+1;
		elseif (n2 > size(newtemp2,1))
			% newtemp2 has been completely scanned
			% The data entirely comes from newtemp1 
			tempedges(i,1:3) = newtemp1(n1,1:3);
			tempedges(i,5:6) = newtemp1(n1,4:5);
			n1 = n1+1;
			i = i+1;
		elseif ((newtemp1(n1,1) == newtemp2(n2,1)) && (newtemp1(n1,2) == newtemp2(n2,2)))
			% the edges from newtemp1 and newtemp2 coincide --> we treat an interior edge
			% The extremities come from newtemp1
			tempedges(i,1:2) = newtemp1(n1,1:2);
			% The index of the cell on the left side (K) of the edge comes from newtemp1
			tempedges(i,3) = newtemp1(n1,3);
			% The index of the cell on the right side (L) of the edge comes from newtemp2
			tempedges(i,4) = newtemp2(n2,3);
			% The indices of preceeding and following vertices in K come from newtemp1
			tempedges(i,5:6) = newtemp1(n1,4:5);
			% The indices of preceeding and following vertices in L come from newtemp1
			tempedges(i,7:8) = newtemp2(n2,4:5);
			i = i+1;
			n1 = n1+1;
			n2 = n2+1;
		elseif ((newtemp1(n1,1) == newtemp2(n2,1)) && (newtemp1(n1,2) < newtemp2(n2,2)))
			% The line from newtemp1 preceeds the line from newtemp2 in the sorting convention
			% --> The data entirely comes from newtemp1
			tempedges(i,1:3) = newtemp1(n1,1:3);
			tempedges(i,5:6) = newtemp1(n1,4:5);
			i = i+1;
			n1 = n1+1;
		elseif (newtemp1(n1,1) < newtemp2(n2,1))
			% The line from newtemp1 preceeds the line from newtemp2 in the sorting convention
			% --> The data entirely comes from newtemp1
			tempedges(i,1:3) = newtemp1(n1,1:3);
			tempedges(i,5:6) = newtemp1(n1,4:5);
			i = i+1;
			n1 = n1+1;
		else
			% The line from newtemp2 preceeds the line from newtemp1 in the sorting convention
			% --> The data entirely comes from newtemp2
			tempedges(i,1:3) = newtemp2(n2,1:3);
			tempedges(i,5:6) = newtemp2(n2,4:5);
			i = i+1;
			n2 = n2+1;
		end
	end

	clear newtemp1;
	clear newtemp2;

	% Extract the real list of edges
	ne = 1;
	while (~((tempedges(ne,1) == 0) && (tempedges(ne,2) == 0)))
		ne = ne+1;
	end
	ne = ne-1;

	etemp = [tempedges(1:ne,1:2)';
		(-2*(tempedges(1:ne,4) < ones(ne,1))+ones(ne,1))';
		tempedges(1:ne,3:4)'];

	clear tempedges;


	% We correct the orientation of the edges in order to fit with the orientation of cells
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	et = zeros(ne,5);
	switch PARAMETERS.PARALLELIZATION.TOOLBOX
		case {'SERIAL'}
			for i=1:ne
				Kstar = etemp(1,i);
				Lstar = etemp(2,i);
				K = etemp(4,i);
				% We verify if the edge (xK*,xL*) is correctly oriented in cell K
				% If it is not the case, we switch the indices K* and L* in the edges list
				edge_is_OK = (((Kstar == t(1,K)) && (Lstar == t(2,K))) || ...
								((Kstar == t(2,K)) && (Lstar == t(3,K))) || ...
								((Kstar == t(3,K)) && (Lstar == t(1,K))));
				et(i,:) = [ ...
								etemp(1,i)*edge_is_OK+etemp(2,i)*(1-edge_is_OK), ...
								etemp(2,i)*edge_is_OK+etemp(1,i)*(1-edge_is_OK), ...
								etemp(3:5,i)'];
			end
		case {'PCT'}
			parfor i=1:ne
				Kstar = etemp(1,i);
				Lstar = etemp(2,i);
				K = etemp(4,i);
				% We verify if the edge (xK*,xL*) is correctly oriented in cell K
				% If it is not the case, we switch the indices K* and L* in the edges list
				edge_is_OK = (((Kstar == t(1,K)) && (Lstar == t(2,K))) || ...
								((Kstar == t(2,K)) && (Lstar == t(3,K))) || ...
								((Kstar == t(3,K)) && (Lstar == t(1,K))));
				et(i,:) = [ ...
								etemp(1,i)*edge_is_OK+etemp(2,i)*(1-edge_is_OK), ...
								etemp(2,i)*edge_is_OK+etemp(1,i)*(1-edge_is_OK), ...
								etemp(3:5,i)'];
			end
		otherwise
			error(sprintf('Unknown value of PARAMETERS.PARALLELIZATION.TOOLBOX\n'));
	end
	e = et';

	clear etemp;

	nbe = 0;
	for i=1:ne
		if (e(5,i) < 1)
			nbe = nbe+1;
		end
	end
	be = zeros(5,nbe);
	k = 1;
	for i=1:ne
		if (e(5,i) < 1)
			be(:,k) = e(:,i);
			be(5,k) = i;
			k = k+1;
		end
	end


	% The full list of edges is built. We now add the description of cells by edges.
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	tenosort = zeros(3,nt);
	countedge_triangle = zeros(1,nt);
	for i=1:ne
		K = e(4,i);
		L = e(5,i);
		countedge_triangle(K) = countedge_triangle(K)+1;
		tenosort(countedge_triangle(K),K) = i;
		if (L > 0)
			countedge_triangle(L) = countedge_triangle(L)+1;
			tenosort(countedge_triangle(L),L) = -i;
		end
	end

	clear countedge_triangle;

	tet = zeros(nt,3);
	switch PARAMETERS.PARALLELIZATION.TOOLBOX
		case {'SERIAL'}
			for ke=1:nt
				for ip=1:3
					inode = t(ip,ke);
					for ie=1:3
						iedge = abs(tenosort(ie,ke));
						if (~(inode == e(1,iedge)) && ~(inode == e(2,iedge)))
							% inode is the vertex that is opposite to the edge iedge in the triangle ke
							tet(ke,ip) = iedge;
						end
					end
				end
			end
		case {'PCT'}
			parfor ke=1:nt
				for ip=1:3
					inode = t(ip,ke);
					for ie=1:3
						iedge = abs(tenosort(ie,ke));
						if (~(inode == e(1,iedge)) && ~(inode == e(2,iedge)))
							% inode is the vertex that is opposite to the edge iedge in the triangle ke
							tet(ke,ip) = iedge;
						end
					end
				end
			end
		otherwise
			error(sprintf('Unknown value of PARAMETERS.PARALLELIZATION.TOOLBOX\n'));
	end
	te = tet';

	clear tenosort;
	clear tet;

	if (min(te) <= 0)
		error(sprintf('Problem while building the descriptor by edges of triangles.\n'));
	end

end