identify_boundsegments_P1B.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[BS_top, BS_bottom, BS_left, BS_right] = identify_boundsegments_P1B(MESH, DOMAIN)
% INPUT
% MESH Mesh metadata
% DOMAIN Domain metadata
%
% OUTPUT
% BS_top List of P1B bound vertices on the top part of the domain
% BS_bottom List of P1B bound vertices on the bottom part of the domain
% BS_left List of P1B bound vertices on the left part of the domain
% BS_right List of P1B bound vertices on the right part of the domain
BS_top = [];
BS_bottom = [];
BS_left = [];
BS_right = [];
cprec = 1e-8;
switch DOMAIN.GEOMETRY
case {'RECTANGLE'}
if isfield(MESH, 'BS_top_p2')
% No need to build BS_* lists with loops
BS_top = [MESH.BS_top_p2(1,:); MESH.BS_top_p2(3,:)];
BS_bottom = [MESH.BS_bottom_p2(1,:); MESH.BS_bottom_p2(3,:)];
BS_left = [MESH.BS_left_p2(1,:); MESH.BS_left_p2(3,:)];
BS_right = [MESH.BS_right_p2(1,:); MESH.BS_right_p2(3,:)];
else
Lx = DOMAIN.XMAX-DOMAIN.XMIN;
Ly = DOMAIN.YMAX-DOMAIN.YMIN;
% Loop on bound edges
for ibe1b=1:MESH.nbe1b
% Find the index of the bound edge in the full list of edges
ie1b = MESH.be1b(5,ibe1b);
% Find the index of the extremity vertices
iv1 = MESH.e1bfull(1,ie1b);
iv2 = MESH.e1bfull(2,ie1b);
% Find the unique triangle K such that sigma is an edge of K
K = MESH.e1bfull(4,ie1b);
if (((iv1 == MESH.t1b(1,K)) && (iv2 == MESH.t1b(2,K))) || ...
((iv2 == MESH.t1b(1,K)) && (iv1 == MESH.t1b(2,K))))
locseg = [MESH.t1b(1,K); MESH.t1b(2,K)];
elseif (((iv1 == MESH.t1b(2,K)) && (iv2 == MESH.t1b(3,K))) || ...
((iv2 == MESH.t1b(2,K)) && (iv1 == MESH.t1b(3,K))))
locseg = [MESH.t1b(2,K); MESH.t1b(3,K)];
elseif (((iv1 == MESH.t1b(3,K)) && (iv2 == MESH.t1b(1,K))) || ...
((iv2 == MESH.t1b(3,K)) && (iv1 == MESH.t1b(1,K))))
locseg = [MESH.t1b(3,K); MESH.t1b(1,K)];
else
error('Problem while identifying the segments on boundary');
end
% Find the coordinates of these vertices
v1(1:2) = MESH.p1b(1:2,iv1);
v2(1:2) = MESH.p1b(1:2,iv2);
if (max(abs(v1(2)-DOMAIN.YMAX)/Ly, abs(v2(2)-DOMAIN.YMAX)/Ly) < cprec)
% The bound edge is located on top bound
BS_top = [BS_top, locseg];
elseif (max(abs(v1(2)-DOMAIN.YMIN)/Ly, abs(v2(2)-DOMAIN.YMIN)/Ly) < cprec)
% The bound edge is located on bottom bound
BS_bottom = [BS_bottom, locseg];
elseif (max(abs(v1(1)-DOMAIN.XMIN)/Lx, abs(v2(1)-DOMAIN.XMIN)/Lx) < cprec)
% The bound edge is located on left bound
BS_left = [BS_left, locseg];
elseif (max(abs(v1(1)-DOMAIN.XMAX)/Lx, abs(v2(1)-DOMAIN.XMAX)/Lx) < cprec)
% The bound edge is located on right bound
BS_right = [BS_right, locseg];
else
% The bound edge cannot be localized on the rectangle boundary
error('Problem while localizing a bound edge\n');
end
end
end
case {'DIHEDRON'}
if isfield(MESH, 'BS_top_p2')
% No need to build BS_* lists with loops
BS_top = [MESH.BS_top_p2(1,:); MESH.BS_top_p2(3,:)];
BS_bottom = [MESH.BS_bottom_p2(1,:); MESH.BS_bottom_p2(3,:)];
BS_left = [MESH.BS_left_p2(1,:); MESH.BS_left_p2(3,:)];
BS_right = [MESH.BS_right_p2(1,:); MESH.BS_right_p2(3,:)];
else
Lx = DOMAIN.XMAX-DOMAIN.XMIN;
Ly = DOMAIN.YMIN-DOMAIN.YMIN;
% Loop on bound edges
for ibe1b=1:MESH.nbe1b
% Find the index of the bound edge in the full list of edges
ie1b = MESH.be1b(5,ibe1b);
% Find the index of the extremity vertices
iv1 = MESH.e1bfull(1,ie1b);
iv2 = MESH.e1bfull(2,ie1b);
% Find the unique triangle K such that sigma is an edge of K
K = MESH.e1bfull(4,ie1b);
if (((iv1 == MESH.t1b(1,K)) && (iv2 == MESH.t1b(2,K))) || ...
((iv2 == MESH.t1b(1,K)) && (iv1 == MESH.t1b(2,K))))
locseg = [MESH.t1b(1,K); MESH.t1b(2,K)];
elseif (((iv1 == MESH.t1b(2,K)) && (iv2 == MESH.t1b(3,K))) || ...
((iv2 == MESH.t1b(2,K)) && (iv1 == MESH.t1b(3,K))))
locseg = [MESH.t1b(2,K); MESH.t1b(3,K)];
elseif (((iv1 == MESH.t1b(3,K)) && (iv2 == MESH.t1b(1,K))) || ...
((iv2 == MESH.t1b(3,K)) && (iv1 == MESH.t1b(1,K))))
locseg = [MESH.t1b(3,K); MESH.t1b(1,K)];
else
error('Problem while identifying the segments on boundary');
end
% Find the coordinates of these vertices
v1(1:2) = MESH.p1b(1:2,iv1);
v2(1:2) = MESH.p1b(1:2,iv2);
a = (DOMAIN.YMIN-DOMAIN.YAXIS_ANGLE1) / (DOMAIN.XAXIS_ANGLE2-DOMAIN.XAXIS_ANGLE1);
b = DOMAIN.YAXIS_ANGLE1-a*DOMAIN.XAXIS_ANGLE1;
test_top = (max(abs(v1(2)-DOMAIN.YMAX)/Ly, abs(v2(2)-DOMAIN.YMAX)/Ly) < cprec);
test_high_plateau = ((max(abs(v1(2)-DOMAIN.YAXIS_ANGLE1)/Ly, ...
abs(v2(2)-DOMAIN.YAXIS_ANGLE1)/Ly) < cprec) && ...
(min(v1(1),v2(1)) >= DOMAIN.XMIN-cprec*Lx) && ...
(max(v1(1),v2(1)) <= DOMAIN.XAXIS_ANGLE1+cprec*Lx));
test_low_plateau = ((max(abs(v1(2)-DOMAIN.YMIN)/Ly, ...
abs(v2(2)-DOMAIN.YMIN)/Ly) < cprec) && ...
(min(v1(1),v2(1)) >= DOMAIN.XAXIS_ANGLE2-cprec*Lx) && ...
(max(v1(1),v2(1)) <= DOMAIN.XMAX+cprec*Lx));
test_fall = ((max(abs(v1(2)-a*v1(1)-b)/Ly, abs(v2(2)-a*v2(1)-b)/Ly) < cprec) && ...
(min(v1(1),v2(1)) >= DOMAIN.XAXIS_ANGLE1-cprec*Lx) && ...
(max(v1(1),v2(1)) <= DOMAIN.XAXIS_ANGLE2+cprec*Lx));
test_bottom = test_low_plateau || test_high_plateau || test_fall;
test_left = (max(abs(v1(1)-DOMAIN.XMIN)/Lx, abs(v2(1)-DOMAIN.XMIN)/Lx) < cprec);
test_right = (max(abs(v1(1)-DOMAIN.XMAX)/Lx, abs(v2(1)-DOMAIN.XMAX)/Lx) < cprec);
if test_top
% The bound edge is located on top bound
BS_top = [BS_top, locseg];
elseif test_bottom
% The bound edge is located on bottom bound
BS_bottom = [BS_bottom, locseg];
elseif test_left
% The bound edge is located on left bound
BS_left = [BS_left, locseg];
elseif test_right
% The bound edge is located on right bound
BS_right = [BS_right, locseg];
else
% The bound edge cannot be localized on the dihedron boundary
error('Problem while localizing a bound edge');
end
end
end
case {'STEP'}
if isfield(MESH, 'BS_top_p2')
% No need to build BS_* lists with loops
BS_top = [MESH.BS_top_p2(1,:); MESH.BS_top_p2(3,:)];
BS_bottom = [MESH.BS_bottom_p2(1,:); MESH.BS_bottom_p2(3,:)];
BS_left = [MESH.BS_left_p2(1,:); MESH.BS_left_p2(3,:)];
BS_right = [MESH.BS_right_p2(1,:); MESH.BS_right_p2(3,:)];
else
Lx = DOMAIN.XMAX-DOMAIN.XMIN;
Ly = DOMAIN.YMAX-DOMAIN.YMIN;
% Loop on bound edges
for ibe1b=1:MESH.nbe1b
% Find the index of the bound edge in the full list of edges
ie1b = MESH.be1b(5,ibe1b);
% Find the index of the extremity vertices
iv1 = MESH.e1bfull(1,ie1b);
iv2 = MESH.e1bfull(2,ie1b);
% Find the unique triangle K such that sigma is an edge of K
K = MESH.e1bfull(4,ie1b);
if (((iv1 == MESH.t1b(1,K)) && (iv2 == MESH.t1b(2,K))) || ...
((iv2 == MESH.t1b(1,K)) && (iv1 == MESH.t1b(2,K))))
locseg = [MESH.t1b(1,K); MESH.t1b(2,K)];
elseif (((iv1 == MESH.t1b(2,K)) && (iv2 == MESH.t1b(3,K))) || ...
((iv2 == MESH.t1b(2,K)) && (iv1 == MESH.t1b(3,K))))
locseg = [MESH.t1b(2,K); MESH.t1b(3,K)];
elseif (((iv1 == MESH.t1b(3,K)) && (iv2 == MESH.t1b(1,K))) || ...
((iv2 == MESH.t1b(3,K)) && (iv1 == MESH.t1b(1,K))))
locseg = [MESH.t1b(3,K); MESH.t1b(1,K)];
else
error('Problem while identifying the segments on boundary');
end
% Find the coordinates of these vertices
v1(1:2) = MESH.p1b(1:2,iv1);
v2(1:2) = MESH.p1b(1:2,iv2);
test_top = (max(abs(v1(2)-DOMAIN.YMAX)/Ly, abs(v2(2)-DOMAIN.YMAX)/Ly) < cprec);
test_left = (max(abs(v1(1)-DOMAIN.XMIN)/Lx, abs(v2(1)-DOMAIN.XMIN)/Lx) < cprec);
test_right = (max(abs(v1(1)-DOMAIN.XMAX)/Lx, abs(v2(1)-DOMAIN.XMAX)/Lx) < cprec);
test_fall = ((max(abs(v1(1)-DOMAIN.XSTEP)/Lx, abs(v2(1)-DOMAIN.XSTEP)/Lx) < cprec) && ...
(min(v1(2),v2(2)) >= DOMAIN.YMIN) && (max(v1(2),v2(2)) <= DOMAIN.YSTEP));
switch DOMAIN.STEP_POSITION
case {'LEFT'}
test_plateau_left = ((min(v1(1),v2(1)) >= DOMAIN.XMIN) && ...
(max(v1(1),v2(1)) <= DOMAIN.XSTEP) && ...
(max(abs(v1(2)-DOMAIN.YSTEP)/Ly, abs(v2(2)-DOMAIN.YSTEP)/Ly) < cprec));
test_plateau_right = ((min(v1(1),v2(1)) >= DOMAIN.XSTEP) && ...
(max(v1(1),v2(1)) <= DOMAIN.XMAX) && ...
(max(abs(v1(2)-DOMAIN.YMIN)/Ly, abs(v2(2)-DOMAIN.YMIN)/Ly) < cprec));
case {'RIGHT'}
test_plateau_left = ((min(v1(1),v2(1)) >= DOMAIN.XMIN) && ...
(max(v1(1),v2(1)) <= DOMAIN.XSTEP) && ...
(max(abs(v1(2)-DOMAIN.YMIN)/Ly, abs(v2(2)-DOMAIN.YMIN)/Ly) < cprec));
test_plateau_right = ((min(v1(1),v2(1)) >= DOMAIN.XSTEP) && ...
(max(v1(1),v2(1)) <= DOMAIN.XMAX) && ...
(max(abs(v1(2)-DOMAIN.YSTEP)/Ly, abs(v2(2)-DOMAIN.YSTEP)/Ly) < cprec));
otherwise
error('Unknown position of step in STEP geometry');
end
test_bottom = test_fall || test_plateau_left || test_plateau_right;
if test_top
% The bound edge is located on top bound
BS_top = [BS_top, locseg];
elseif test_bottom
% The bound edge is located on bottom bound
BS_bottom = [BS_bottom, locseg];
elseif test_left
% The bound edge is located on left bound
BS_left = [BS_left, locseg];
elseif test_right
% The bound edge is located on right bound
BS_right = [BS_right, locseg];
else
% The bound edge cannot be localized on the dihedron boundary
error('Problem while localizing a bound edge');
end
end
end
case {'CIRCLE', 'FROM_FILE'}
% Nothing to do
otherwise
error('Unknown domain geometry');
end
end
