mesh_stars_fv_part1.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 [normal,volume,volume_P2,perimeter,volume_barycenter,barycenter_patchEF] = mesh_stars_fv_part1(MESH)
% INPUT
% MESH An incomplete mesh metadata
%
% OUTPUT
% normal Normal unit vector to the edges
% volume Control Volume surface for density with P1
% volume_P2 Control Volume surface for density with P2
% perimeter Control Volume perimeter for density with P1
% volume_barycenter Control Volume barycenter coordonates
% barycenter_patchEF Barycenter computation with FE patch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% VOLUME and NORMAL quantities %%%%%%
volume = zeros(1,MESH.np1);
normal = zeros(6,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))];
for j = 1:3
volume (MESH.t1(j,ke)) = volume(MESH.t1(j,ke)) + abs(det(De))/6;
normal (2*j-1,ke) = -(MESH.p2(2,MESH.t2(j+3,ke)) - MESH.p2(2,MESH.t2(j,ke)))/3;
normal (2*j ,ke) = (MESH.p2(1,MESH.t2(j+3,ke)) - MESH.p2(1,MESH.t2(j,ke)))/3;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% PERIMETER quantity of each cell %%%%%%
perimeter = zeros(1,MESH.np1);
for ke = 1:MESH.nt1
for j = 1:3
xe = (MESH.p2(1,MESH.t2(j+3,ke)) - MESH.p2(1,MESH.t2(j,ke)))/3;
ye = (MESH.p2(2,MESH.t2(j+3,ke)) - MESH.p2(2,MESH.t2(j,ke)))/3;
nor = sqrt(xe^2+ye^2);
jp1 = mod(j ,3)+1;
jm1 = mod(jp1,3)+1;
neup = MESH.t1(jp1,ke);
neum = MESH.t1(jm1,ke);
perimeter(neup) = perimeter(neup) + nor;
perimeter(neum) = perimeter(neum) + nor;
end
end
% For a node on the boundary, add the bondary part of the perimeter
for ia = 1:MESH.ne1full % edges number
if (MESH.e1full(3,ia) ~= 0)
xe = (MESH.p1(1,MESH.e1full(1,ia)) - MESH.p1(1,MESH.e1full(2,ia)))/2;
ye = (MESH.p1(2,MESH.e1full(1,ia)) - MESH.p1(2,MESH.e1full(2,ia)))/2;
nor = sqrt(xe^2+ye^2);
perimeter(MESH.e1full(1,ia)) = perimeter(MESH.e1full(1,ia)) + nor;
perimeter(MESH.e1full(2,ia)) = perimeter(MESH.e1full(2,ia)) + nor;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Volume_P2 : Control volume surface for P2 density %%%%%%
% PTSOM = sub-triangles nodes
ptsom(1,1) = 1;
ptsom(1,2) = 6;
ptsom(1,3) = 5;
ptsom(2,1) = 6;
ptsom(2,2) = 2;
ptsom(2,3) = 4;
ptsom(3,1) = 5;
ptsom(3,2) = 4;
ptsom(3,3) = 3;
ptsom(4,1) = 4;
ptsom(4,2) = 5;
ptsom(4,3) = 6;
volume_P2 = zeros(1,MESH.np2);
for ke=1:MESH.nt1 % loop on elements
for stri=1:4 % loop on sub triangles
% sub-triangle nodes
n1 = ptsom(stri,1);
n2 = ptsom(stri,2);
n3 = ptsom(stri,3);
% elemental Jacobian = each sub-triangle surface
De = [MESH.p2(2,MESH.t2(n3,ke))-MESH.p2(2,MESH.t2(n1,ke)), ...
MESH.p2(2,MESH.t2(n1,ke))-MESH.p2(2,MESH.t2(n2,ke)); ...
MESH.p2(1,MESH.t2(n1,ke))-MESH.p2(1,MESH.t2(n3,ke)), ...
MESH.p2(1,MESH.t2(n2,ke))-MESH.p2(1,MESH.t2(n1,ke))];
for j = 1:3
jsom = ptsom(stri,j);
volume_P2(MESH.t2(jsom,ke)) = volume_P2(MESH.t2(jsom,ke)) + abs(det(De))/6;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% barycenter computation with Control Volumes %%%%
volume_barycenter = zeros(2,MESH.np1);
bari_tri1 = zeros(6,MESH.nt1);
bari_tri2 = zeros(6,MESH.nt1);
for ke = 1:MESH.nt1
% cell barycenter
bari_elem(1) = (MESH.p1(1,MESH.t1(1,ke))+MESH.p1(1,MESH.t1(2,ke))+MESH.p1(1,MESH.t1(3,ke)))/3;
bari_elem(2) = (MESH.p1(2,MESH.t1(1,ke))+MESH.p1(2,MESH.t1(2,ke))+MESH.p1(2,MESH.t1(3,ke)))/3;
for j = 1:3
jp1 = mod(j ,3)+1+3;
jm1 = mod(jp1,3)+1+3;
% surface of 2 KE sub-triangles built by linking the midle of
% an edge to the barycenter
De_tr1 = [bari_elem(2)-MESH.p2(2,MESH.t2(j,ke)), ...
MESH.p2(2,MESH.t2(j,ke))-MESH.p2(2,MESH.t2(jp1,ke)); ...
MESH.p2(1,MESH.t2(j,ke))-bari_elem(1), ...
MESH.p2(1,MESH.t2(jp1,ke))-MESH.p2(1,MESH.t2(j,ke))];
De_tr2 = [bari_elem(2)-MESH.p2(2,MESH.t2(j,ke)), ...
MESH.p2(2,MESH.t2(j,ke))-MESH.p2(2,MESH.t2(jm1,ke)); ...
MESH.p2(1,MESH.t2(j,ke))-bari_elem(1), ...
MESH.p2(1,MESH.t2(jm1,ke))-MESH.p2(1,MESH.t2(j,ke))];
A1 = abs(det(De_tr1))/2;
A2 = abs(det(De_tr2))/2;
bari_tri1(2*j-1,ke) = (MESH.p2(1,MESH.t2(j,ke)) + MESH.p2(1,MESH.t2(jp1,ke)) + bari_elem(1))/3;
bari_tri1(2*j ,ke) = (MESH.p2(2,MESH.t2(j,ke)) + MESH.p2(2,MESH.t2(jp1,ke)) + bari_elem(2))/3;
bari_tri2(2*j-1,ke) = (MESH.p2(1,MESH.t2(j,ke)) + MESH.p2(1,MESH.t2(jm1,ke)) + bari_elem(1))/3;
bari_tri2(2*j ,ke) = (MESH.p2(2,MESH.t2(j,ke)) + MESH.p2(2,MESH.t2(jm1,ke)) + bari_elem(2))/3;
volume_barycenter(1,MESH.t1(j,ke)) = volume_barycenter(1,MESH.t1(j,ke)) + ...
bari_tri1(2*j-1,ke)*A1 + bari_tri2(2*j-1,ke)*A2;
volume_barycenter(2,MESH.t1(j,ke)) = volume_barycenter(2,MESH.t1(j,ke)) + ...
bari_tri1(2*j,ke)*A1 + bari_tri2(2*j,ke)*A2;
end
end
volume_barycenter(1,:) = volume_barycenter(1,:)./volume;
volume_barycenter(2,:) = volume_barycenter(2,:)./volume;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% barycenter computation with FE patch %%%%
barycenter_patchEF = zeros(2,MESH.np1);
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))];
% barycenter of each triangle
bari_elem(1) = (MESH.p1(1,MESH.t1(1,ke))+MESH.p1(1,MESH.t1(2,ke))+MESH.p1(1,MESH.t1(3,ke)))/3;
bari_elem(2) = (MESH.p1(2,MESH.t1(1,ke))+MESH.p1(2,MESH.t1(2,ke))+MESH.p1(2,MESH.t1(3,ke)))/3;
for j = 1:3
barycenter_patchEF(1,MESH.t1(j,ke)) = barycenter_patchEF(1,MESH.t1(j,ke)) + bari_elem(1)*abs(det(De))/2;
barycenter_patchEF(2,MESH.t1(j,ke)) = barycenter_patchEF(2,MESH.t1(j,ke)) + bari_elem(2)*abs(det(De))/2;
end
end
barycenter_patchEF(1,:) = barycenter_patchEF(1,:)./(volume*3);
barycenter_patchEF(2,:) = barycenter_patchEF(2,:)./(volume*3);
end
