edge_descriptor.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
