modify_adjacent_cells_border.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 [adjacent_cells,flag_shadow_cells,corners] = modify_adjacent_cells_border(MESH, PARAMETERS)
% INPUT
% MESH An incomplete mesh metadata
% PARAMETERS List of simulation parameters
%
% OUTPUT
% adjacent_cells Contains the upstream and downstream cells of each edge
% flag_shadow_cells Flag each shadow upstream or downstream triangle
% corners List of four corners of a rectangular domain
%% Compute the list of boundary P1-nodes and identify each boundary
% This function locates the upstream and downstream cells of each edge,
% only for triangles which share a boundary node
% for RECTANGLE or CIRCLE domains
Nboundary = MESH.be1(1,:);
Idboundary = MESH.be1(5,:);
nb_border_node = length(Nboundary);
%% Compute the list of four corners of a rectangular domain
corners = [];
if strcmp(PARAMETERS.DOMAIN.GEOMETRY,'RECTANGLE')
for i=1:nb_border_node
if ((MESH.p1(1,Nboundary(i)) == PARAMETERS.DOMAIN.XMIN && ...
MESH.p1(2,Nboundary(i)) == PARAMETERS.DOMAIN.YMIN) || ...
(MESH.p1(1,Nboundary(i)) == PARAMETERS.DOMAIN.XMAX && ...
MESH.p1(2,Nboundary(i)) == PARAMETERS.DOMAIN.YMIN) || ...
(MESH.p1(1,Nboundary(i)) == PARAMETERS.DOMAIN.XMIN && ...
MESH.p1(2,Nboundary(i)) == PARAMETERS.DOMAIN.YMAX) || ...
(MESH.p1(1,Nboundary(i)) == PARAMETERS.DOMAIN.XMAX && ...
MESH.p1(2,Nboundary(i)) == PARAMETERS.DOMAIN.YMAX))
corners = [corners,Nboundary(i)];
end
end
end
%% Compute the number of triangles which share a boundary node
triangles_border_node = zeros(1,nb_border_node);
for ke = 1:MESH.nt1
for j=1:3
js = MESH.t1(j,ke);
[member,pos] = ismember(js,Nboundary);
if member
triangles_border_node(pos) = triangles_border_node(pos)+1;
end
end
end
max_triangles_border_node = max(triangles_border_node);
%% Take note of the number of triangles which share a boundary node
num_triangles_border_node = zeros(max_triangles_border_node,nb_border_node);
for ke = 1:MESH.nt1
for j=1:3
js = MESH.t1(j,ke);
[member,pos] = ismember(js,Nboundary);
ind = 1;
if member
while (num_triangles_border_node(ind,pos) > 0)
ind = ind+1;
end
num_triangles_border_node(ind,pos) = ke;
end
end
end
%% Modify in adjacent_cells the upstream and downstream cells for each edge of each triangle which shares a boundary node
mask_triangle = zeros(1,MESH.nt1);
flag_shadow_cells = zeros(MESH.nt1,3,2);
adjacent_cells = MESH.adjacent_cells;
switch PARAMETERS.DOMAIN.GEOMETRY
case {'RECTANGLE','DIHEDRON'}
% When a triangle ke contains a boundary node, we consider the symetric of the vector (ddx,ddy)
% w.r.t. the vertical or horizontal boundary and we flag the "shadow" triangle
% Rmq: there not exists a shadow triangle for the four corners of the domain
for c=1:nb_border_node
for r=1:triangles_border_node(c)
ke = num_triangles_border_node(r,c); % global number of triangle
if ~mask_triangle(ke) % this is a triangle not yet visited
for j=1:3
% Compute the middle point of the segment [xG,xsigma] where xG is the gravity
% center of the triangle and xsigma is the middle point of sigma
% Name it Cj
xCj = (sum(MESH.p2(1,MESH.t2(1:3,ke)))/3. + MESH.p2(1,MESH.t2(j+3,ke)))/2.;
yCj = (sum(MESH.p2(2,MESH.t2(1:3,ke)))/3. + MESH.p2(2,MESH.t2(j+3,ke)))/2.;
% Identify the vertices that are the extremities of the edge
% (counter-clockwise orientation)
% Name them Aj+ and Aj-
jp1 = mod(j,3)+1;
jm1 = mod(jp1,3)+1;
Ajp = MESH.t1(jp1,ke);
Ajm = MESH.t1(jm1,ke);
% Verify if Aj+ is a boundary point
[member,pos] = ismember(Ajp,Nboundary);
if member
% Compute the vector AjpCj and normalize it
AjpCjx = xCj-MESH.p1(1,Ajp);
AjpCjy = yCj-MESH.p1(2,Ajp);
nor = sqrt(AjpCjx^2+AjpCjy^2);
AjpCjx = AjpCjx/nor;
AjpCjy = AjpCjy/nor;
pdscalmin = 2;
ind_triangle = 0;
for it=1:triangles_border_node(pos)
nloc = 1;
nb_triangle = num_triangles_border_node(it,pos);
while (MESH.t1(nloc,nb_triangle) ~= Ajp)
nloc = nloc+1;
end
node_face1 = MESH.t2(nloc+3,nb_triangle);
xv = MESH.p2(1,node_face1)-MESH.p1(1,Ajp); % the renumbering is not the same for p1 and p2
yv = MESH.p2(2,node_face1)-MESH.p1(2,Ajp);
nor = sqrt(xv^2+yv^2);
xv = xv/nor;
yv = yv/nor;
pdscal = AjpCjx*xv+AjpCjy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
shadow = 0;
if ((Idboundary(pos) == 2) || (Idboundary(pos) == 4)) % Ajp in vertical boundary
pdscal = -AjpCjx*xv+AjpCjy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
shadow = 1;
elseif ((Idboundary(pos) == 1) || (Idboundary(pos) == 3)) % Ajp in horizontal boundary
pdscal = AjpCjx*xv-AjpCjy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
shadow = 2;
end
end
adjacent_cells(ke,j,1) = ind_triangle;
flag_shadow_cells(ke,j,1) = shadow; % flag the shadow triangle ke
end
% Verify if Aj- is a boundary point
[member,pos] = ismember(Ajm,Nboundary);
if member
% Compute the vector AjmCj and normalize it
AjmCjx = xCj-MESH.p1(1,Ajm);
AjmCjy = yCj-MESH.p1(2,Ajm);
nor = sqrt(AjmCjx^2+AjmCjy^2);
AjmCjx = AjmCjx/nor;
AjmCjy = AjmCjy/nor;
pdscalmin = 2;
ind_triangle = 0;
for it=1:triangles_border_node(pos)
nloc = 1;
nb_triangle = num_triangles_border_node(it,pos);
while (MESH.t1(nloc,nb_triangle) ~= Ajm)
nloc = nloc+1;
end
node_face2 = MESH.t2(nloc+3,nb_triangle);
xv = MESH.p2(1,node_face2)-MESH.p1(1,Ajm); % the renumbering is not the same for p1 and p2
yv = MESH.p2(2,node_face2)-MESH.p1(2,Ajm);
nor = sqrt(xv^2+yv^2);
xv = xv/nor;
yv = yv/nor;
pdscal = AjmCjx*xv+AjmCjy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
shadow = 0;
if ((Idboundary(pos) == 2) || (Idboundary(pos) == 4)) % Ajm in vertical boundary
pdscal = -AjmCjx*xv+AjmCjy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
shadow = 1;
elseif ((Idboundary(pos) == 1) || (Idboundary(pos) == 3)) % Ajm in horizontal boundary
pdscal = AjmCjx*xv-AjmCjy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
shadow = 2;
end
end
adjacent_cells(ke,j,2) = ind_triangle;
flag_shadow_cells(ke,j,2) = shadow; % flag the shadow triangle ke
end
end
mask_triangle(ke) = 1;
end
end
end
case{'CIRCLE'}
% When a triangle ke contains a boundary node, we consider the symmetric of
% the vector (ddx,ddy) w.r.t. the tangent to the circle and we flag the "shadow" triangle
for c=1:nb_border_node
for r=1:triangles_border_node(c)
ke = num_triangles_border_node(r,c); % global number of triangle
if ~mask_triangle(ke) % this is a triangle not yet visited
for j=1:3
% Compute the middle point of the segment [xG,xsigma] where xG is the gravity
% center of the triangle and xsigma is the middle point of sigma
% Name it Cj
xCj = (sum(MESH.p2(1,MESH.t2(1:3,ke)))/3 + MESH.p2(1,MESH.t2(j+3,ke)))/2.;
yCj = (sum(MESH.p2(2,MESH.t2(1:3,ke)))/3 + MESH.p2(2,MESH.t2(j+3,ke)))/2.;
% Identify the vertices that are the extremities of the edge
% (counter-clockwise orientation)
% Name them Aj+ and Aj-
jp1 = mod(j,3)+1;
jm1 = mod(jp1,3)+1;
Ajp = MESH.t1(jp1,ke);
Ajm = MESH.t1(jm1,ke);
% Verify if Aj+ is a boundary point
[member,pos] = ismember(Ajp,Nboundary);
if member
% Compute the vector AjpCj and normalize it
AjpCjx = xCj-MESH.p1(1,Ajp);
AjpCjy = yCj-MESH.p1(2,Ajp);
nor = sqrt(AjpCjx^2+ddy^2);
AjpCjx = AjpCjx/nor;
AjpCjy = AjpCjy/nor;
xb = MESH.p1(1,Ajp)-PARAMETERS.DOMAIN.XCENTER;
yb = MESH.p1(2,Ajp)-PARAMETERS.DOMAIN.YCENTER;
nor = sqrt(xb^2+yb^2);
cos_a = xb/nor;
sin_a = yb/nor;
pdscalmin = 2;
ind_triangle = 0;
for it=1:triangles_border_node(pos)
nloc = 1;
nb_triangle = num_triangles_border_node(it,pos);
while (MESH.t1(nloc,nb_triangle) ~=neup)
nloc = nloc+1;
end
node_face1 = MESHt2(nloc+3,nb_triangle);
xv = MESH.p2(1,node_face1)-MESH.p1(1,Ajp); % the renumbering is not the same for p1 and p2
yv = MESH.p2(2,node_face1)-MESH.p1(2,Ajp);
nor = sqrt(xv^2+yv^2);
xv = xv/nor;
yv = yv/nor;
pdscal = AjpCjx*xv+AjpCjy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
shadow = 0;
sym_AjpCjx = (sin_a^2-cos_a^2)*AjpCjx-2.*sin_a*cos_a*AjpCjy; % symmetric of the vector AjpCj w.r.t. the tangent to the circle
sym_AjpCjy = (cos_a^2-sin_a^2)*AjpCjy-2.*sin_a*cos_a*AjpCjx;
pdscal = sym_AjpCjx*xv+sym_AjpCjy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
shadow = 1;
end
adjacent_cells(ke,j,1) = ind_triangle;
flag_shadow_cells(ke,j,1) = shadow; % flag the shadow triangle ke
end
% Verify if Aj- is a boundary point
[member,pos] = ismember(Ajm,Nboundary);
if member
AjmCjx = xCj - MESH.p1(1,Ajm);
AjmCjy = yCj - MESH.p1(2,Ajm);
nor = sqrt(AjmCjx^2+AjmCjy^2);
AjmCjx = AjmCjx/nor;
AjmCjy = AjmCjy/nor;
xb = MESH.p1(1,Ajm)-PARAMETERS.DOMAIN.XCENTER;
yb = MESH.p1(2,Ajm)-PARAMETERS.DOMAIN.YCENTER;
nor = sqrt(xb^2+yb^2);
cos_a = xb/nor;
sin_a = yb/nor;
pdscalmin = 2;
ind_triangle = 0;
for it=1:triangles_border_node(pos)
nloc = 1;
nb_triangle = num_triangles_border_node(it,pos);
while (MESH.t1(nloc,nb_triangle) ~= Ajm)
nloc = nloc+1;
end
node_face2 = MESH.t2(nloc+3,nb_triangle);
xv = MESH.p2(1,node_face2)-MESH.p1(1,Ajm); % the renumbering is not the same for p1 and p2
yv = MESH.p2(2,node_face2)-MESH.p1(2,Ajm);
nor = sqrt(xv^2+yv^2);
xv = xv/nor;
yv = yv/nor;
pdscal = AjmCjx*xv+AjmCjy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
shadow = 0;
sym_AjmCjx = (sin_a^2-cos_a^2)*AjmCjx-2.*sin_a*cos_a*AjmCjy; % symmetric of the vector AjmCj w.r.t. the tangent to the circle
sym_AjmCjy = (cos_a^2-sin_a^2)*AjmCjy-2.*sin_a*cos_a*AjmCjx;
pdscal = sym_AjmCjx*xv+sym_AjmCjy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
shadow = 1;
end
adjacent_cells(ke,j,2) = ind_triangle;
flag_shadow_cells(ke,j,2) = shadow; % flag the shadow triangle ke
end
end
mask_triangle(ke) = 1;
end
end
end
case {'FROM_FILE'}
error(sprintf('The routine for computing shadow cells does not handle geometries provided in external files\n'));
otherwise
error(sprintf('Unknown domain geometry.\n'));
end
end
