structmesh_dihedron_2d.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[p,t] = structmesh_dihedron_2d(PARAMETERS)
% INPUT
% PARAMETERS List of simulation parameters
%
% OUTPUT
% p P1 nodes
% t Triangles described by P1 nodes
% Mesh a rectangle first
%%%%%%%%%%%%%%%%%%%%%%%%
NBSEG_X = PARAMETERS.MESH.NBSEG_PERUNIT_X * (PARAMETERS.DOMAIN.XMAX-PARAMETERS.DOMAIN.XMIN);
NBSEG_Y = PARAMETERS.MESH.NBNODES_INBL_Y-1 + PARAMETERS.MESH.NBSEG_OUTBL_Y;
np = (NBSEG_X+1)*(NBSEG_Y+1);
nt = 2*NBSEG_X*NBSEG_Y;
ptemp = zeros(np,4);
ttemp = zeros(nt,5);
x = linspace(PARAMETERS.DOMAIN.XMIN, PARAMETERS.DOMAIN.XMAX, NBSEG_X+1);
yinbl = linspace(PARAMETERS.DOMAIN.YAXIS_ANGLE1, ...
PARAMETERS.DOMAIN.YAXIS_ANGLE1+PARAMETERS.MESH.HEIGHT_BL, PARAMETERS.MESH.NBNODES_INBL_Y);
dyout = (PARAMETERS.DOMAIN.YMAX-PARAMETERS.DOMAIN.YAXIS_ANGLE1-PARAMETERS.MESH.HEIGHT_BL) / ...
PARAMETERS.MESH.NBSEG_OUTBL_Y;
youtbl = linspace(PARAMETERS.DOMAIN.YAXIS_ANGLE1 + PARAMETERS.MESH.HEIGHT_BL + dyout, ...
PARAMETERS.DOMAIN.YMAX, PARAMETERS.MESH.NBSEG_OUTBL_Y);
y = [yinbl, youtbl];
% Coordinates of nodes
for i=1:NBSEG_X+1
for j=1:NBSEG_Y+1
num = (j-1)*(NBSEG_X+1) + i;
ptemp(num,1) = num;
ptemp(num,2) = x(i);
ptemp(num,3) = y(j);
if ((i == 1) || (i == NBSEG_X+1) || (j == 1) || (j == NBSEG_Y+1))
ptemp(num,4) = 1;
end
end
end
% Triangles
switch PARAMETERS.MESH.TRIANGLES_ORIENTATION
case {'DIAGONAL'}
ke = 1;
for i=1:NBSEG_X
for j=1:NBSEG_Y
ttemp(ke,1) = ke;
ttemp(ke,2) = (j-1)*(NBSEG_X+1) + i;
ttemp(ke,3) = (j-1)*(NBSEG_X+1) + i+1;
ttemp(ke,4) = j*(NBSEG_X+1) + i+1;
ttemp(ke,5) = 0;
ke = ke+1;
ttemp(ke,1) = ke;
ttemp(ke,2) = (j-1)*(NBSEG_X+1) + i;
ttemp(ke,3) = j*(NBSEG_X+1) + i+1;
ttemp(ke,4) = j*(NBSEG_X+1) + i;
ttemp(ke,5) = 0;
ke = ke+1;
end
end
case {'CROSS'}
if (mod(NBSEG_X,2) ~= 0)
error('The number of edges on top and bottom bounds must be a multiple of 2.\n');
end
ke = 1;
for j=1:floor(NBSEG_Y/2)
for i=1:floor(NBSEG_X/2)
% bottom left quadrant, diagonal : bottom left - upper right
is = (2*j-2) * (NBSEG_X+1) + 2*i-1;
ttemp(ke,1) = ke;
ttemp(ke,2) = is+1;
ttemp(ke,3) = is+NBSEG_X+2;
ttemp(ke,4) = is;
ttemp(ke,5) = 0;
ke = ke+1;
ttemp(ke,1) = ke;
ttemp(ke,2) = is+NBSEG_X+1;
ttemp(ke,3) = is;
ttemp(ke,4) = is+NBSEG_X+2;
ttemp(ke,5) = 0;
ke = ke+1;
% bottom right quadrant, diagonal : upper left - bottom right
is = (2*j-2) * (NBSEG_X+1) + 2*i;
ttemp(ke,1) = ke;
ttemp(ke,2) = is;
ttemp(ke,3) = is+1;
ttemp(ke,4) = is+NBSEG_X+1;
ttemp(ke,5) = 0;
ke = ke+1;
ttemp(ke,1) = ke;
ttemp(ke,2) = is+NBSEG_X+2;
ttemp(ke,3) = is+NBSEG_X+1;
ttemp(ke,4) = is+1;
ttemp(ke,5) = 0;
ke = ke+1;
% upper left quadrant, diagonal : upper left - bottom right
is = (2*j-1) * (NBSEG_X+1) + 2*i-1;
ttemp(ke,1) = ke;
ttemp(ke,2) = is;
ttemp(ke,3) = is+1;
ttemp(ke,4) = is+NBSEG_X+1;
ttemp(ke,5) = 0;
ke = ke+1;
ttemp(ke,1) = ke;
ttemp(ke,2) = is+NBSEG_X+2;
ttemp(ke,3) = is+NBSEG_X+1;
ttemp(ke,4) = is+1;
ttemp(ke,5) = 0;
ke = ke+1;
% upper right quadrant, diagonal : bottom left - upper right
is = (2*j-1) * (NBSEG_X+1) + 2*i;
ttemp(ke,1) = ke;
ttemp(ke,2) = is+1;
ttemp(ke,3) = is+NBSEG_X+2;
ttemp(ke,4) = is;
ttemp(ke,5) = 0;
ke = ke+1;
ttemp(ke,1) = ke;
ttemp(ke,2) = is+NBSEG_X+1;
ttemp(ke,3) = is;
ttemp(ke,4) = is+NBSEG_X+2;
ttemp(ke,5) = 0;
ke = ke+1;
end
end
if (mod(NBSEG_Y,2) ~= 0) % Add another layer of triangles
j=ceil(NBSEG_Y/2);
for i=1:floor(NBSEG_X/2)
% bottom left quadrant, diagonal : bottom left - upper right
is = (2*j-2) * (NBSEG_X+1) + 2*i-1;
ttemp(ke,1) = ke;
ttemp(ke,2) = is+1;
ttemp(ke,3) = is+NBSEG_X+2;
ttemp(ke,4) = is;
ttemp(ke,5) = 0;
ke = ke+1;
ttemp(ke,1) = ke;
ttemp(ke,2) = is+NBSEG_X+1;
ttemp(ke,3) = is;
ttemp(ke,4) = is+NBSEG_X+2;
ttemp(ke,5) = 0;
ke = ke+1;
% bottom right quadrant, diagonal : upper left - bottom right
is = (2*j-2) * (NBSEG_X+1) + 2*i;
ttemp(ke,1) = ke;
ttemp(ke,2) = is;
ttemp(ke,3) = is+1;
ttemp(ke,4) = is+NBSEG_X+1;
ttemp(ke,5) = 0;
ke = ke+1;
ttemp(ke,1) = ke;
ttemp(ke,2) = is+NBSEG_X+2;
ttemp(ke,3) = is+NBSEG_X+1;
ttemp(ke,4) = is+1;
ttemp(ke,5) = 0;
ke = ke+1;
end
end
otherwise
error('Wrong triangle orientation.\n')
end
% Modify y-component of nodes in the boundary layer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:NBSEG_X+1
for j=1:PARAMETERS.MESH.NBNODES_INBL_Y
ip = (j-1)*(NBSEG_X+1) + i;
ptemp(ip,3) = PARAMETERS.DOMAIN.YAXIS_ANGLE1 + ((ptemp(ip,3)-PARAMETERS.DOMAIN.YAXIS_ANGLE1) / ...
PARAMETERS.MESH.HEIGHT_BL)^PARAMETERS.MESH.PROG_GEOM_BL_Y * PARAMETERS.MESH.HEIGHT_BL;
end
end
% Transform into a dihedron
% Case of boundary layer
%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:NBSEG_X+1
if (ptemp(i,2) <= PARAMETERS.DOMAIN.XAXIS_ANGLE1)
ybottom = PARAMETERS.DOMAIN.YAXIS_ANGLE1;
elseif (ptemp(i,2) >= PARAMETERS.DOMAIN.XAXIS_ANGLE2)
ybottom = PARAMETERS.DOMAIN.YMIN;
else
ybottom = PARAMETERS.DOMAIN.YAXIS_ANGLE1 + (ptemp(i,2)-PARAMETERS.DOMAIN.XAXIS_ANGLE1) / ...
(PARAMETERS.DOMAIN.XAXIS_ANGLE2-PARAMETERS.DOMAIN.XAXIS_ANGLE1) * ...
(PARAMETERS.DOMAIN.YMIN-PARAMETERS.DOMAIN.YAXIS_ANGLE1);
end
for j=1:PARAMETERS.MESH.NBNODES_INBL_Y
ip = (j-1)*(NBSEG_X+1)+i;
ptemp(ip,3) = ybottom + ptemp(ip,3) - PARAMETERS.DOMAIN.YAXIS_ANGLE1;
end
end
% Case out of boundary layer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=PARAMETERS.MESH.NBNODES_INBL_Y+1:NBSEG_Y+1
yup = ptemp((j-1)*(NBSEG_X+1)+1,3);
ydown = PARAMETERS.DOMAIN.YMIN + PARAMETERS.MESH.HEIGHT_BL + (j-PARAMETERS.MESH.NBNODES_INBL_Y) * ...
(PARAMETERS.DOMAIN.YMAX-PARAMETERS.DOMAIN.YMIN-PARAMETERS.MESH.HEIGHT_BL) / ...
PARAMETERS.MESH.NBSEG_OUTBL_Y;
for i=((PARAMETERS.DOMAIN.XAXIS_ANGLE1-PARAMETERS.DOMAIN.XMIN)*PARAMETERS.MESH.NBSEG_PERUNIT_X+1):(NBSEG_X+1)
ip = (j-1)*(NBSEG_X+1) + i;
if (ptemp(ip,2) >= PARAMETERS.DOMAIN.XAXIS_ANGLE2)
ptemp(ip,3) = ydown;
else
ptemp(ip,3) = yup + (ptemp(ip,2)-PARAMETERS.DOMAIN.XAXIS_ANGLE1) / ...
(PARAMETERS.DOMAIN.XAXIS_ANGLE2-PARAMETERS.DOMAIN.XAXIS_ANGLE1)*(ydown-yup);
end
end
end
p = ptemp(:,2:3)';
t = ttemp(:,2:4)';
end
