mesh_squares_fv.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[horizontal_edges, vertical_edges, volume, perimeter] = mesh_squares_fv(MESH, PARAMETERS)
% INPUT
% MESH An incomplete mesh metadata
% PARAMETERS List of simulation parameters
%
% OUTPUT
% horizontal_edges Contains the i,j index of the horizontal edge of each triangle
% vertical_edges Contains the i,j index of the verical edge of each triangle
% volume Control Volume surface for density with P1
% perimeter Control Volume barycenter coordonates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Construction of HORIZONTAL_EDGES and VERTICAL_EDGES arrays %%%%
%%%% similar to connectivity array building in rect_2d.m %%%%
horizontal_edges = zeros(2,MESH.nt1);
vertical_edges = zeros(2,MESH.nt1);
ie = 0;
switch PARAMETERS.MESH.TRIANGLES_ORIENTATION
case {'DIAGONAL'}
isx = PARAMETERS.MESH.NBSEG_X/2;
switch PARAMETERS.MODEL
case {'NSDV'}
switch PARAMETERS.TESTCASE
case {'RTIN'}
jsyb = 3*PARAMETERS.MESH.NBSEG_Y/4;
case {'DROP'}
jsyb = MESH.NBSEG_Y/4;
case {'CAEN', 'EXAC', 'EXAC2', 'EXACTNEU'}
jsyb = PARAMETERS.MESH.NBSEG_Y/2;
otherwise
error(sprintf('Unknown test case.\n'));
end
otherwise
error(sprintf('Unknown model.\n'));
end
% bottom left quadrant, diagonal : bottom left - upper right
for j = 1:jsyb
for i = 1:isx
ie = ie + 1;
horizontal_edges(1,ie) = i;
horizontal_edges(2,ie) = j;
vertical_edges(1,ie) = i+1;
vertical_edges(2,ie) = j;
ie = ie + 1;
horizontal_edges(1,ie) = i;
horizontal_edges(2,ie) = j+1;
vertical_edges(1,ie) = i;
vertical_edges(2,ie) = j;
end
end
% bottom right quadrant, diagonal : upper left - bottom right
for j = 1:jsyb
for i = isx+1:PARAMETERS.MESH.NBSEG_X;
ie = ie + 1;
horizontal_edges(1,ie) = i;
horizontal_edges(2,ie) = j;
vertical_edges(1,ie) = i;
vertical_edges(2,ie) = j;
ie = ie + 1;
horizontal_edges(1,ie) = i;
horizontal_edges(2,ie) = j+1;
vertical_edges(1,ie) = i+1;
vertical_edges(2,ie) = j;
end
end
% upper left quadrant, diagonal : upper left - bottom right
for j = jsyb+1:PARAMETERS.MESH.NBSEG_Y
for i = 1:isx
ie = ie + 1;
horizontal_edges(1,ie) = i;
horizontal_edges(2,ie) = j;
vertical_edges(1,ie) = i;
vertical_edges(2,ie) = j;
ie = ie + 1;
horizontal_edges(1,ie) = i;
horizontal_edges(2,ie) = j+1;
vertical_edges(1,ie) = i+1;
vertical_edges(2,ie) = j;
end
end
% upper right quadrant, diagonal : bottom left - upper right
for j = jsyb+1:PARAMETERS.MESH.NBSEG_Y
for i = isx+1:PARAMETERS.MESH.NBSEG_X
ie = ie + 1;
horizontal_edges(1,ie) = i;
horizontal_edges(2,ie) = j;
vertical_edges(1,ie) = i+1;
vertical_edges(2,ie) = j;
ie = ie + 1;
horizontal_edges(1,ie) = i;
horizontal_edges(2,ie) = j+1;
vertical_edges(1,ie) = i;
vertical_edges(2,ie) = j;
end
end
case {'CROSS'}
for j = 1:PARAMETERS.MESH.NBSEG_Y/2
for i = 1:PARAMETERS.MESH.NBSEG_X/2
% bottom left quadrant, diagonal : bottom left - upper right
ie = ie + 1;
horizontal_edges(1,ie) = 2*i-1;
horizontal_edges(2,ie) = 2*j-1;
vertical_edges(1,ie) = 2*i;
vertical_edges(2,ie) = 2*j-1;
ie = ie + 1;
horizontal_edges(1,ie) = 2*i-1;
horizontal_edges(2,ie) = 2*j;
vertical_edges(1,ie) = 2*i-1;
vertical_edges(2,ie) = 2*j-1;
% bottom right quadrant, diagonal : upper left - bottom right
ie = ie + 1;
horizontal_edges(1,ie) = 2*i;
horizontal_edges(2,ie) = 2*j-1;
vertical_edges(1,ie) = 2*i;
vertical_edges(2,ie) = 2*j-1;
ie = ie + 1;
horizontal_edges(1,ie) = 2*i;
horizontal_edges(2,ie) = 2*j;
vertical_edges(1,ie) = 2*i+1;
vertical_edges(2,ie) = 2*j-1;
% upper left quadrant, diagonal : upper left - bottom right
ie = ie + 1;
horizontal_edges(1,ie) = 2*i-1;
horizontal_edges(2,ie) = 2*j;
vertical_edges(1,ie) = 2*i-1;
vertical_edges(2,ie) = 2*j;
ie = ie + 1;
horizontal_edges(1,ie) = 2*i-1;
horizontal_edges(2,ie) = 2*j+1;
vertical_edges(1,ie) = 2*i;
vertical_edges(2,ie) = 2*j;
% upper right quadrant, diagonal : bottom left - upper right
ie = ie + 1;
horizontal_edges(1,ie) = 2*i;
horizontal_edges(2,ie) = 2*j;
vertical_edges(1,ie) = 2*i+1;
vertical_edges(2,ie) = 2*j;
ie = ie + 1;
horizontal_edges(1,ie) = 2*i;
horizontal_edges(2,ie) = 2*j+1;
vertical_edges(1,ie) = 2*i;
vertical_edges(2,ie) = 2*j;
end
end
otherwise
error(sprintf('Unknown triangle orientation.\n'));
end
Dx = (PARAMETERS.DOMAIN.XMAX-PARAMETERS.DOMAIN.XMIN)/PARAMETERS.MESH.NBSEG_X;
Dy = (PARAMETERS.DOMAIN.YMAX-PARAMETERS.DOMAIN.YMIN)/PARAMETERS.MESH.NBSEG_Y;
volume = (Dx*Dy)*ones(1,MESH.np1);
perimeter = 2.*(Dx+Dy)*ones(1,MESH.np1);
% border volumes are half volumes
for J=2:PARAMETERS.MESH.NBSEG_Y
K1 = (PARAMETERS.MESH.NBSEG_X+1)*(J-1)+1;
K2 = (PARAMETERS.MESH.NBSEG_X+1)*J;
volume(K1) = volume(K1)/2.;
volume(K2) = volume(K2)/2.;
perimeter(K1) = perimeter(K1)-Dx;
perimeter(K2) = perimeter(K2)-Dx;
end
for I=2:PARAMETERS.MESH.NBSEG_X
K1 = I;
K2 = PARAMETERS.MESH.NBSEG_Y*(PARAMETERS.MESH.NBSEG_X+1)+I;
volume(K1) = volume(K1)/2.;
volume(K2) = volume(K2)/2.;
perimeter(K1) = perimeter(K1)-Dy;
perimeter(K2) = perimeter(K2)-Dy;
end
% corner volumes are quarter volumes
volume(1) = Dx*Dy/4;
volume(PARAMETERS.MESH.NBSEG_X+1) = Dx*Dy/4;
volume(PARAMETERS.MESH.NBSEG_Y*(PARAMETERS.MESH.NBSEG_X+1)+1) = Dx*Dy/4;
volume((PARAMETERS.MESH.NBSEG_Y+1)*(PARAMETERS.MESH.NBSEG_X+1)) = Dx*Dy/4;
perimeter(1) = Dx+Dy;
perimeter(PARAMETERS.MESH.NBSEG_X+1) = Dx+Dy;
perimeter(PARAMETERS.MESH.NBSEG_Y*(PARAMETERS.MESH.NBSEG_X+1)+1) = Dx+Dy;
perimeter((PARAMETERS.MESH.NBSEG_Y+1)*(PARAMETERS.MESH.NBSEG_X+1)) = Dx+Dy;
end
