structmesh_rectangle_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_rectangle_2d(PARAMETERS)
% INPUT
% PARAMETERS List of simulation parameters
%
% OUTPUT
% p P1 nodes
% t Triangles described by P1 nodes
nt = PARAMETERS.MESH.NBSEG_X * PARAMETERS.MESH.NBSEG_Y * 2;
np = (PARAMETERS.MESH.NBSEG_X+1) * (PARAMETERS.MESH.NBSEG_Y+1);
t = zeros(3,nt);
p = zeros(2,np);
switch PARAMETERS.MESH.TRIANGLES_ORIENTATION
case {'DIAGONAL'}
if (mod(PARAMETERS.MESH.NBSEG_X,2) ~= 0)
error('The number of edges on North and South bounds must be a multiple of 2.\n');
end
if ((mod(PARAMETERS.MESH.NBSEG_Y,4) ~= 0) && (max(strcmp(PARAMETERS.TESTCASE, {'RTIN','DROP'})) == 1))
error('The number of edges on West and East bounds must be a multiple of 4.\n');
end
isx = PARAMETERS.MESH.NBSEG_X/2;
switch PARAMETERS.MODEL
case {'NSDV'}
switch PARAMETERS.TESTCASE
case {'RTIN'}
jsyb = 3*PARAMETERS.MESH.NBSEG_Y/4;
jsyh = PARAMETERS.MESH.NBSEG_Y/4;
case {'DROP'}
jsyb = PARAMETERS.MESH.NBSEG_Y/4;
jsyh = 3*PARAMETERS.MESH.NBSEG_Y/4;
case {'CAEN', 'EXAC', 'EXAC2', 'EXACNEU'}
jsyb = PARAMETERS.MESH.NBSEG_Y/2;
jsyh = PARAMETERS.MESH.NBSEG_Y/2;
otherwise
error('Unknown test case.\n');
end
case {'NS'}
switch PARAMETERS.TESTCASE
case {'CAEN', 'EXAC', 'EXACNEU', 'POIS', 'NONA', 'GTPSI'}
jsyb = PARAMETERS.MESH.NBSEG_Y/2;
jsyh = PARAMETERS.MESH.NBSEG_Y/2;
otherwise
error('Unknown test case.\n');
end
otherwise
error('Unknown model.\n');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% mesh connectivities array building %%%%
ie = 0;
is0 = 0;
% bottom left quadrant, diagonal : bottom left - upper right
for j = 1:jsyb
for i = 1:isx
is = is0 + (j-1) * (PARAMETERS.MESH.NBSEG_X+1) + i;
ie = ie + 1;
t(1,ie) = is + 1;
t(2,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
t(3,ie) = is;
ie = ie + 1;
t(1,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
t(2,ie) = is;
t(3,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
end
end
is0 = isx;
% bottom right quadrant, diagonal : upper left - bottom right
for j = 1:jsyb
for i = 1:isx
is = is0 + (j-1) * (PARAMETERS.MESH.NBSEG_X+1) + i;
ie = ie + 1;
t(1,ie) = is;
t(2,ie) = is + 1;
t(3,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
ie = ie + 1;
t(1,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
t(2,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
t(3,ie) = is + 1;
end
end
is0 = jsyb * (PARAMETERS.MESH.NBSEG_X+1);
% upper left quadrant, diagonal : upper left - bottom right
for j = 1:jsyh
for i = 1:isx
is = is0 + (j-1) * (PARAMETERS.MESH.NBSEG_X+1) + i;
ie = ie + 1;
t(1,ie) = is;
t(2,ie) = is + 1;
t(3,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
ie = ie + 1;
t(1,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
t(2,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
t(3,ie) = is + 1;
end
end
is0 = jsyb * (PARAMETERS.MESH.NBSEG_X+1) + isx;
% upper right quadrant, diagonal : bottom left - upper right
for j = 1:jsyh
for i = 1:isx
is = is0 + (j-1) * (PARAMETERS.MESH.NBSEG_X+1) + i;
ie = ie + 1;
t(1,ie) = is + 1;
t(2,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
t(3,ie) = is;
ie = ie + 1;
t(1,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
t(2,ie) = is;
t(3,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
end
end
case {'CROSS'}
if (mod(PARAMETERS.MESH.NBSEG_X,2) ~= 0)
error('The number of edges on North and South bounds must be a multiple of 2.\n');
end
if (mod(PARAMETERS.MESH.NBSEG_Y,2) ~= 0)
error('The number of edges on West and East bounds must be a multiple of 2.\n');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% mesh connectivities array building %%%%
ie = 0;
for j = 1:PARAMETERS.MESH.NBSEG_Y/2
for i = 1:PARAMETERS.MESH.NBSEG_X/2
% bottom left quadrant, diagonal : bottom left - upper right
is = (2*j-2) * (PARAMETERS.MESH.NBSEG_X+1) + 2*i-1;
ie = ie + 1;
t(1,ie) = is + 1;
t(2,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
t(3,ie) = is;
ie = ie + 1;
t(1,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
t(2,ie) = is;
t(3,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
% bottom right quadrant, diagonal : upper left - bottom right
is = (2*j-2) * (PARAMETERS.MESH.NBSEG_X+1) + 2*i;
ie = ie + 1;
t(1,ie) = is;
t(2,ie) = is + 1;
t(3,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
ie = ie + 1;
t(1,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
t(2,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
t(3,ie) = is + 1;
% upper left quadrant, diagonal : upper left - bottom right
is = (2*j-1) * (PARAMETERS.MESH.NBSEG_X+1) + 2*i-1;
ie = ie + 1;
t(1,ie) = is;
t(2,ie) = is + 1;
t(3,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
ie = ie + 1;
t(1,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
t(2,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
t(3,ie) = is + 1;
% upper right quadrant, diagonal : bottom left - upper right
is = (2*j-1) * (PARAMETERS.MESH.NBSEG_X+1) + 2*i;
ie = ie + 1;
t(1,ie) = is + 1;
t(2,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
t(3,ie) = is;
ie = ie + 1;
t(1,ie) = is + PARAMETERS.MESH.NBSEG_X+1;
t(2,ie) = is;
t(3,ie) = is + PARAMETERS.MESH.NBSEG_X+2;
end
end
otherwise
error('Wrong triangle orientation.\n')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Nodes coordonates array building %%%%
step_ix = (PARAMETERS.DOMAIN.XMAX - PARAMETERS.DOMAIN.XMIN) / PARAMETERS.MESH.NBSEG_X;
step_iy = (PARAMETERS.DOMAIN.YMAX - PARAMETERS.DOMAIN.YMIN) / PARAMETERS.MESH.NBSEG_Y;
is = 0;
for j = 1:PARAMETERS.MESH.NBSEG_Y+1
for i = 1:PARAMETERS.MESH.NBSEG_X+1
is = is + 1;
p(1,is) = PARAMETERS.DOMAIN.XMIN + step_ix * (i-1);
p(2,is) = PARAMETERS.DOMAIN.YMIN + step_iy * (j-1);
end
end
end
