mesh_stars_fv_part3.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[upstream_cell,downstream_cell,adjacent_cells,ciarlet_ratio] = mesh_stars_fv_part3(MESH)
% INPUT
% MESH An incomplete mesh metadata
%
% OUTPUT
% upstream_cell Locates the downstream cells of each edge
% downstream_cell Locates the downstream cells of each edge
% adjacent_cells An other way to describe upstream and downstream cells of each edge
% ciarlet_ratio Ciarlet_ratio number as defined in Calgaro-Chane-Kane-Creuse-Goudon JCP2010
%% Compute the number of triangles which share a fixed node
% This function locates the upstream and downstream cells of each edge
% An other way to describe upstream and downstream cells of each edge
% Compute the ciarlet_ratio number
triangles_node = zeros(1,MESH.np1);
for ke = 1:MESH.nt1
for j=1:3
js = MESH.t1(j,ke);
triangles_node(js) = triangles_node(js)+1;
end
end
max_triangles_node = max(triangles_node);
%% Take note of the number of triangles which share a fixed node
num_triangles_node = zeros(max_triangles_node,MESH.np1);
for ke = 1:MESH.nt1
for j=1:3
js = MESH.t1(j,ke);
ind = 1;
while (num_triangles_node(ind,js) > 0)
ind = ind+1;
end
num_triangles_node(ind,js) = ke;
end
end
%% Compute the upstream_cell and downstream_cell of each edge
upstream_cell = zeros(1,MESH.ne1full);
downstream_cell = zeros(1,MESH.ne1full);
for ie=1:MESH.ne1full
node1 = MESH.e1full(1,ie);
node2 = MESH.e1full(2,ie);
% (xe,ye) is the normal vector of the edge
xe = MESH.p1(1,node2)-MESH.p1(1,node1);
ye = MESH.p1(2,node2)-MESH.p1(2,node1);
nor = sqrt(xe^2+ye^2);
xe = xe/nor;
ye = ye/nor;
pdscalmin = 2;
ind_triangle = 0;
for it=1:triangles_node(node1)
nloc = 1;
nb_triangle = num_triangles_node(it,node1);
while (MESH.t1(nloc,nb_triangle) ~=node1)
nloc = nloc+1;
end
node_face1 = MESH.t2(nloc+3,nb_triangle);
xv = MESH.p2(1,node_face1)-MESH.p1(1,node1); % the renumbering is not the same for p1 and p2
yv = MESH.p2(2,node_face1)-MESH.p1(2,node1);
nor = sqrt(xv^2+yv^2);
xv = xv/nor;
yv = yv/nor;
pdscal = xe*xv+ye*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
end
upstream_cell(1,ie) = ind_triangle;
pdscalmax = -2;
ind_triangle = 0;
for it=1:triangles_node(node2)
nloc = 1;
nb_triangle = num_triangles_node(it,node2);
while (MESH.t1(nloc,nb_triangle) ~=node2)
nloc = nloc+1;
end
node_face2 = MESH.t2(nloc+3,nb_triangle);
xv = MESH.p2(1,node_face2)-MESH.p1(1,node2); % the renumbering is not the same for p1 and p2
yv = MESH.p2(2,node_face2)-MESH.p1(2,node2);
nor = sqrt(xv^2+yv^2);
xv = xv/nor;
yv = yv/nor;
pdscal = xe*xv+ye*yv;
if (pdscal > pdscalmax)
pdscalmax = pdscal;
ind_triangle = nb_triangle;
end
end
downstream_cell(1,ie) = ind_triangle;
end
%% A second way in order to compute the upstream and downstream cells for each edge of each triangle
adjacent_cells = zeros(MESH.nt1,3,2);
for ke=1:MESH.nt1
for j=1:3
xgj = (sum(MESH.p2(1,MESH.t2(1:3,ke)))/3 + MESH.p2(1,MESH.t2(j+3,ke)))/2;
ygj = (sum(MESH.p2(2,MESH.t2(1:3,ke)))/3 + MESH.p2(2,MESH.t2(j+3,ke)))/2;
jp1 = mod(j ,3)+1;
jm1 = mod(jp1,3)+1;
neup = MESH.t1(jp1,ke);
neum = MESH.t1(jm1,ke);
ddx = xgj-MESH.p1(1,neup);
ddy = ygj-MESH.p1(2,neup);
nor = sqrt(ddx^2+ddy^2);
ddx = ddx/nor;
ddy = ddy/nor;
pdscalmin = 2;
ind_triangle = 0;
for it=1:triangles_node(neup)
nloc = 1;
nb_triangle = num_triangles_node(it,neup);
while (MESH.t1(nloc,nb_triangle) ~= neup)
nloc = nloc+1;
end
node_face1 = MESH.t2(nloc+3,nb_triangle);
xv = MESH.p2(1,node_face1)-MESH.p1(1,neup); % the renumbering is not the same for p1 and p2
yv = MESH.p2(2,node_face1)-MESH.p1(2,neup);
nor = sqrt(xv^2+yv^2);
xv = xv/nor;
yv = yv/nor;
pdscal = ddx*xv+ddy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
end
adjacent_cells(ke,j,1) = ind_triangle;
ddx = xgj-MESH.p1(1,neum);
ddy = ygj-MESH.p1(2,neum);
nor = sqrt(ddx^2+ddy^2);
ddx = ddx/nor;
ddy = ddy/nor;
pdscalmin = 2;
ind_triangle = 0;
for it=1:triangles_node(neum)
nloc = 1;
nb_triangle = num_triangles_node(it,neum);
while (MESH.t1(nloc,nb_triangle) ~= neum)
nloc = nloc+1;
end
node_face2 = MESH.t2(nloc+3,nb_triangle);
xv = MESH.p2(1,node_face2)-MESH.p1(1,neum); % the renumbering is not the same for p1 and p2
yv = MESH.p2(2,node_face2)-MESH.p1(2,neum);
nor = sqrt(xv^2+yv^2);
xv = xv/nor;
yv = yv/nor;
pdscal = ddx*xv+ddy*yv;
if (pdscal < pdscalmin)
pdscalmin = pdscal;
ind_triangle = nb_triangle;
end
end
adjacent_cells(ke,j,2) = ind_triangle;
end
end
%% Compute the ciarlet_ratio number
% ciarlet_ratio = max_{j \in triangles} max_{k \in triangles sharing a node with j} r_j/h_k
% r_j being the radius of the circombscribed circle in triangle j
% h_k being the smallest height of the triangle k
radius = zeros(1,MESH.nt1);
min_height = zeros(1,MESH.nt1);
for ke=1:MESH.nt1
De = [MESH.p1(2,MESH.t1(3,ke))-MESH.p1(2,MESH.t1(1,ke)), ...
MESH.p1(2,MESH.t1(1,ke))-MESH.p1(2,MESH.t1(2,ke)); ...
MESH.p1(1,MESH.t1(1,ke))-MESH.p1(1,MESH.t1(3,ke)), ...
MESH.p1(1,MESH.t1(2,ke))-MESH.p1(1,MESH.t1(1,ke))];
Je = abs(det(De));
radius(ke) = 1.0/(2*Je);
min_height(ke) = 1.e+10;
for j=1:3
ia = abs(MESH.t1e(j,ke));
radius(ke) = radius(ke)*MESH.edges_length(ia);
min_height(ke) = min(min_height(ke), Je/MESH.edges_length(ia));
end
end
ciarlet_ratio=0.0;
ratio_rj_hk=0.0;
for ke=1:MESH.nt1
for j=1:3
js = MESH.t1(j,ke);
for it=1:triangles_node(js)
nb_triangle = num_triangles_node(it,js);
ratio_rj_hk = max(ratio_rj_hk, radius(ke) / min_height(nb_triangle));
end
end
ciarlet_ratio = max(ciarlet_ratio,ratio_rj_hk);
end
end
