build_geometry.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[geom] = build_geometry(PARAMETERS)
% INPUT
% PARAMETERS List of simulation parameters
%
% OUTPUT
% geom Domain geometry descriptor for Matlab PDE Toolbox
switch PARAMETERS.DOMAIN.GEOMETRY
case {'RECTANGLE'}
% Define the vertices of the rectangle
v1 = [PARAMETERS.DOMAIN.XMIN; PARAMETERS.DOMAIN.YMIN];
v2 = [PARAMETERS.DOMAIN.XMAX; PARAMETERS.DOMAIN.YMIN];
v3 = [PARAMETERS.DOMAIN.XMAX; PARAMETERS.DOMAIN.YMAX];
v4 = [PARAMETERS.DOMAIN.XMIN; PARAMETERS.DOMAIN.YMAX];
% Domain index at the left of each segment, i.e. the domain index (fixed to 1)
domain_id_left = 1;
% Domain index at the right of each segment, i.e. the exterior (so it is 0)
domain_id_right = 0;
% Define the lines
line1 = [2; v1(1); v2(1); v1(2); v2(2); domain_id_left; domain_id_right];
line2 = [2; v2(1); v3(1); v2(2); v3(2); domain_id_left; domain_id_right];
line3 = [2; v3(1); v4(1); v3(2); v4(2); domain_id_left; domain_id_right];
line4 = [2; v4(1); v1(1); v4(2); v1(2); domain_id_left; domain_id_right];
% Define the geometry
geom = [line1, line2, line3, line4];
case {'CIRCLE'}
% Define the center
c = [PARAMETERS.DOMAIN.XCENTER; PARAMETERS.DOMAIN.YCENTER];
% Consider shift points for defining the perimeter
v1 = c + [PARAMETERS.DOMAIN.RADIUS; 0];
v2 = c + [0; PARAMETERS.DOMAIN.RADIUS];
v3 = c - [PARAMETERS.DOMAIN.RADIUS; 0];
v4 = c - [0; PARAMETERS.DOMAIN.RADIUS];
% Domain index at the left of the perimeter segment, i.e. the domain index (fixed to 1)
domain_id_left = 1;
% Domain index at the right of the perimeter segment, i.e. the exterior (so it is 0)
domain_id_right = 0;
% Define the arcs
arc1 = [1; v1(1); v2(1); v1(2); v2(2); domain_id_left; domain_id_right; c; PARAMETERS.DOMAIN.RADIUS];
arc2 = [1; v2(1); v3(1); v2(2); v3(2); domain_id_left; domain_id_right; c; PARAMETERS.DOMAIN.RADIUS];
arc3 = [1; v3(1); v4(1); v3(2); v4(2); domain_id_left; domain_id_right; c; PARAMETERS.DOMAIN.RADIUS];
arc4 = [1; v4(1); v1(1); v4(2); v1(2); domain_id_left; domain_id_right; c; PARAMETERS.DOMAIN.RADIUS];
% Define the geometry
geom = [arc1, arc2, arc3, arc4];
case {'DIHEDRON'}
v1 = [PARAMETERS.DOMAIN.XMIN ; PARAMETERS.DOMAIN.YAXIS_ANGLE1];
v2 = [PARAMETERS.DOMAIN.XAXIS_ANGLE1; PARAMETERS.DOMAIN.YAXIS_ANGLE1];
v3 = [PARAMETERS.DOMAIN.XAXIS_ANGLE2; PARAMETERS.DOMAIN.YMIN];
v4 = [PARAMETERS.DOMAIN.XMAX ; PARAMETERS.DOMAIN.YMIN];
v5 = [PARAMETERS.DOMAIN.XMAX ; PARAMETERS.DOMAIN.YMAX];
v6 = [PARAMETERS.DOMAIN.XMIN ; PARAMETERS.DOMAIN.YMAX];
if isfield(PARAMETERS.PHYSICAL, 'INJECTION_XLIM')
if (numel(PARAMETERS.PHYSICAL.INJECTION_XLIM) == 2)
a = (PARAMETERS.DOMAIN.YMIN-PARAMETERS.DOMAIN.YAXIS_ANGLE1) / ...
(PARAMETERS.DOMAIN.XAXIS_ANGLE2-PARAMETERS.DOMAIN.XAXIS_ANGLE1);
v7 = [PARAMETERS.PHYSICAL.INJECTION_XLIM{1}, ...
a*(PARAMETERS.PHYSICAL.INJECTION_XLIM{1}-PARAMETERS.DOMAIN.XAXIS_ANGLE1) ...
+ PARAMETERS.DOMAIN.YAXIS_ANGLE1];
v8 = [PARAMETERS.PHYSICAL.INJECTION_XLIM{2}, ...
a*(PARAMETERS.PHYSICAL.INJECTION_XLIM{2}-PARAMETERS.DOMAIN.XAXIS_ANGLE1) ...
+ PARAMETERS.DOMAIN.YAXIS_ANGLE1];
end
end
% Domain index at the left of each segment, i.e. the domain index (fixed to 1)
domain_id_left = 1;
% Domain index at the right of each segment, i.e. the exterior (so it is 0)
domain_id_right = 0;
if isfield(PARAMETERS.PHYSICAL, 'INJECTION_XLIM')
if (numel(PARAMETERS.PHYSICAL.INJECTION_XLIM) == 2)
% Define the lines
line1 = [2; v1(1); v2(1); v1(2); v2(2); domain_id_left; domain_id_right];
line2 = [2; v2(1); v7(1); v2(2); v7(2); domain_id_left; domain_id_right];
line3 = [2; v7(1); v8(1); v7(2); v8(2); domain_id_left; domain_id_right];
line4 = [2; v8(1); v3(1); v8(2); v3(2); domain_id_left; domain_id_right];
line5 = [2; v3(1); v4(1); v3(2); v4(2); domain_id_left; domain_id_right];
line6 = [2; v4(1); v5(1); v4(2); v5(2); domain_id_left; domain_id_right];
line7 = [2; v5(1); v6(1); v5(2); v6(2); domain_id_left; domain_id_right];
line8 = [2; v6(1); v1(1); v6(2); v1(2); domain_id_left; domain_id_right];
% Define the geometry
geom = [line1, line2, line3, line4, line5, line6, line7, line8];
else
% Define the lines
line1 = [2; v1(1); v2(1); v1(2); v2(2); domain_id_left; domain_id_right];
line2 = [2; v2(1); v3(1); v2(2); v3(2); domain_id_left; domain_id_right];
line3 = [2; v3(1); v4(1); v3(2); v4(2); domain_id_left; domain_id_right];
line4 = [2; v4(1); v5(1); v4(2); v5(2); domain_id_left; domain_id_right];
line5 = [2; v5(1); v6(1); v5(2); v6(2); domain_id_left; domain_id_right];
line6 = [2; v6(1); v1(1); v6(2); v1(2); domain_id_left; domain_id_right];
% Define the geometry
geom = [line1, line2, line3, line4, line5, line6];
end
else
% Define the lines
line1 = [2; v1(1); v2(1); v1(2); v2(2); domain_id_left; domain_id_right];
line2 = [2; v2(1); v3(1); v2(2); v3(2); domain_id_left; domain_id_right];
line3 = [2; v3(1); v4(1); v3(2); v4(2); domain_id_left; domain_id_right];
line4 = [2; v4(1); v5(1); v4(2); v5(2); domain_id_left; domain_id_right];
line5 = [2; v5(1); v6(1); v5(2); v6(2); domain_id_left; domain_id_right];
line6 = [2; v6(1); v1(1); v6(2); v1(2); domain_id_left; domain_id_right];
% Define the geometry
geom = [line1, line2, line3, line4, line5, line6];
end
case {'STEP'}
switch PARAMETERS.DOMAIN.STEP_POSITION
case {'LEFT'}
v1 = [PARAMETERS.DOMAIN.XMIN; PARAMETERS.DOMAIN.YSTEP];
v2 = [PARAMETERS.DOMAIN.XSTEP; PARAMETERS.DOMAIN.YSTEP];
v3 = [PARAMETERS.DOMAIN.XSTEP; PARAMETERS.DOMAIN.YMIN];
v4 = [PARAMETERS.DOMAIN.XMAX; PARAMETERS.DOMAIN.YMIN];
v5 = [PARAMETERS.DOMAIN.XMAX; PARAMETERS.DOMAIN.YMAX];
v6 = [PARAMETERS.DOMAIN.XMIN; PARAMETERS.DOMAIN.YMAX];
case {'RIGHT'}
v1 = [PARAMETERS.DOMAIN.XMIN; PARAMETERS.DOMAIN.YMIN];
v2 = [PARAMETERS.DOMAIN.XSTEP; PARAMETERS.DOMAIN.YMIN];
v3 = [PARAMETERS.DOMAIN.XSTEP; PARAMETERS.DOMAIN.YSTEP];
v4 = [PARAMETERS.DOMAIN.XMAX; PARAMETERS.DOMAIN.YSTEP];
v5 = [PARAMETERS.DOMAIN.XMAX; PARAMETERS.DOMAIN.YMAX];
v6 = [PARAMETERS.DOMAIN.XMIN; PARAMETERS.DOMAIN.YMAX];
otherwise
error('Unknown position of step in STEP domain geometry');
end
% Domain index at the left of each segment, i.e. the domain index (fixed to 1)
domain_id_left = 1;
% Domain index at the right of each segment, i.e. the exterior (so it is 0)
domain_id_right = 0;
% Define the lines
line1 = [2; v1(1); v2(1); v1(2); v2(2); domain_id_left; domain_id_right];
line2 = [2; v2(1); v3(1); v2(2); v3(2); domain_id_left; domain_id_right];
line3 = [2; v3(1); v4(1); v3(2); v4(2); domain_id_left; domain_id_right];
line4 = [2; v4(1); v5(1); v4(2); v5(2); domain_id_left; domain_id_right];
line5 = [2; v5(1); v6(1); v5(2); v6(2); domain_id_left; domain_id_right];
line6 = [2; v6(1); v1(1); v6(2); v1(2); domain_id_left; domain_id_right];
% Define the geometry
geom = [line1, line2, line3, line4, line5, line6];
otherwise
error('Unknown geometry: only RECTANGLE, CIRCLE, DIHEDRON or STEP geometries are handled.\n');
end
end
