build_mesh.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[MESH, RENUMBERING, DOMAIN] = build_mesh(PARAMETERS, oldp, oldt)
% INPUT
% PARAMETERS List of simulation parameters
% oldp List of P1 nodes coming from an external source (optional)
% oldt List of P1 triangles coming from an external source (optional)
%
% OUTPUTS
% MESH An almost complete mesh metadata
% RENUMBERING Renumbering permutation matrices
% DOMAIN An almost complete domain metadata
RENUMBERING = '';
MESH = '';
DOMAIN = '';
if (exist('oldp') && exist('oldt'))
useoldpt = true;
elseif (~exist('oldp') && ~exist('oldt'))
useoldpt = false;
else
error('Problem with call to build_mesh');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1 - Generate the list of nodes, triangles and edges %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Loading mesh... ');
if ~useoldpt
[MESH.p1, MESH.t1, MESH.hmax, MESH.hmin] = load_mesh(PARAMETERS);
else
MESH.p1 = oldp;
MESH.t1 = oldt;
% Compute hmax and hmin as the maximal and minimal length of an edge in the mesh
MESH.hmax = 0.;
MESH.hmin = max(max(MESH.p1))-min(min(MESH.p1));
for i=1:size(MESH.t1,2)
v1 = MESH.p1(:,MESH.t1(1,i));
v2 = MESH.p1(:,MESH.t1(2,i));
v3 = MESH.p1(:,MESH.t1(3,i));
le1 = sqrt((v1(1)-v2(1))^2+(v1(2)-v2(2))^2);
le2 = sqrt((v2(1)-v3(1))^2+(v2(2)-v3(2))^2);
le3 = sqrt((v3(1)-v1(1))^2+(v3(2)-v1(2))^2);
MESH.hmax = max([MESH.hmax,le1,le2,le3]);
MESH.hmin = min([MESH.hmin,le1,le2,le3]);
end
end
fprintf('done\n');
fprintf('Extracting the edges... ');
[MESH.e1full, MESH.be1, MESH.t1e] = edge_descriptor(MESH.t1, PARAMETERS);
MESH.e2full = MESH.e1full;
MESH.e1bfull = MESH.e1full;
MESH.be2 = MESH.be1;
MESH.be1b = MESH.be1;
MESH.t2e = MESH.t1e;
MESH.t1be = MESH.t1e;
fprintf('done\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2 - Add the number of nodes and the number of triangles %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MESH.nt1 = size(MESH.t1,2);
MESH.np1 = size(MESH.p1,2);
MESH.ne1full = size(MESH.e1full,2);
MESH.nbe1 = size(MESH.be1,2);
MESH.ne2full = size(MESH.e2full,2);
MESH.nbe2 = size(MESH.be2,2);
MESH.ne1bfull = size(MESH.e1bfull,2);
MESH.nbe1b = size(MESH.be1b,2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3- If necessary, complete the definition of the domain geometry and rescale the mesh %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ((strcmp(PARAMETERS.MESH.GENERATION, 'FROM_FILE') ...
&& ~strcmp(PARAMETERS.DOMAIN.GEOMETRY, 'FROM_FILE')) || useoldpt)
fprintf('Extracting the domain geometry... ');
DOMAIN = extract_domain_geometry_from_mesh(MESH, PARAMETERS);
fprintf('done\n');
if (~useoldpt && strcmp(PARAMETERS.MESH.RESCALE, 'YES'))
fprintf('Rescaling the mesh... ');
[MESH.p1, MESH.hmax, MESH.hmin, DOMAIN] = rescale_mesh(MESH, DOMAIN, PARAMETERS);
fprintf('done\n');
end
else
DOMAIN = PARAMETERS.DOMAIN;
%PHYSICAL = PARAMETERS.PHYSICAL;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4 - Add P1B and P2 nodes and triangles %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Computing Finite Element mesh metadata... ');
[MESH.p2, MESH.np2, MESH.t2, MESH.nt2] = add_mdmesh_FE_P2(MESH);
[MESH.p1b, MESH.np1b, MESH.t1b, MESH.nt1b] = add_mdmesh_FE_P1B(MESH);
fprintf('done\n');
%%%%%%%%%%%%%%%%%%%
% 5 - Renumbering %
%%%%%%%%%%%%%%%%%%%
if strcmp(PARAMETERS.MESH.RENUMBERING, 'NONE')
disp('No renumbering in the mesh');
RENUMBERING.perm_p1 = 1:MESH.np1;
RENUMBERING.perm_p1b = 1:MESH.np1b;
RENUMBERING.perm_p2 = 1:MESH.np2;
RENUMBERING.pinv_p1 = 1:MESH.np1;
RENUMBERING.pinv_p1b = 1:MESH.np1b;
RENUMBERING.pinv_p2 = 1:MESH.np2;
else
fprintf('Renumbering in the mesh... ');
% Permutation of the vertices obtained using P1, P1B and P2 DOF
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RENUMBERING = renumbering_P1P1BP2(MESH, PARAMETERS);
RENUMBERING.pinv_p1 = zeros(size(RENUMBERING.perm_p1 ));
RENUMBERING.pinv_p1b = zeros(size(RENUMBERING.perm_p1b));
RENUMBERING.pinv_p2 = zeros(size(RENUMBERING.perm_p2 ));
RENUMBERING.pinv_p1 (RENUMBERING.perm_p1 ) = 1:MESH.np1;
RENUMBERING.pinv_p1b(RENUMBERING.perm_p1b) = 1:MESH.np1b;
RENUMBERING.pinv_p2 (RENUMBERING.perm_p2 ) = 1:MESH.np2;
% Construction of new arrays using the permutations of the vertices
% No permutation for the triangle's number or the edge's number
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MESH.p1 (:,RENUMBERING.perm_p1 ) = MESH.p1;
MESH.p2 (:,RENUMBERING.perm_p2 ) = MESH.p2;
MESH.p1b(:,RENUMBERING.perm_p1b) = MESH.p1b;
for ke=1:MESH.nt1
MESH.t1 (1:3,ke) = RENUMBERING.perm_p1(MESH.t1(1:3,ke));
MESH.t2 (1:6,ke) = RENUMBERING.perm_p2(MESH.t2(1:6,ke));
MESH.t1b(1:4,ke) = RENUMBERING.perm_p1b(MESH.t1b(1:4,ke));
end
for ie=1:MESH.nbe1
MESH.be1(1:2,ie) = RENUMBERING.perm_p1(MESH.be1(1:2,ie));
MESH.be2(1:2,ie) = RENUMBERING.perm_p2(MESH.be2(1:2,ie));
MESH.be1b(1:2,ie) = RENUMBERING.perm_p1b(MESH.be1b(1:2,ie));
end
for ie =1:MESH.ne1full
MESH.e1full(1:2,ie) = RENUMBERING.perm_p1(MESH.e1full(1:2,ie));
MESH.e2full(1:2,ie) = RENUMBERING.perm_p2(MESH.e2full(1:2,ie));
MESH.e1bfull(1:2,ie) = RENUMBERING.perm_p1b(MESH.e1bfull(1:2,ie));
end
fprintf('done\n');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 6- Add Finite element specific descriptors %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Computing Finite Element subtriangles... ');
MESH.st1b = add_subtri_FE_P1B(MESH);
MESH.st2 = add_subtri_FE_P2(MESH);
fprintf('done\n');
fprintf('Identifying the boundary nodes... ');
[MESH.BS_top_p2, MESH.BS_bottom_p2, MESH.BS_left_p2, MESH.BS_right_p2] = ...
identify_boundsegments_P2(MESH, DOMAIN);
[MESH.BS_top_p1b, MESH.BS_bottom_p1b, MESH.BS_left_p1b, MESH.BS_right_p1b] = ...
identify_boundsegments_P1B(MESH, DOMAIN);
fprintf('done\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 7 - Add Finite volume specific descriptors %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Computing Finite Volume mesh metadata... ');
if isfield(PARAMETERS,'FV')
% Add metadata on triangles described by edges
[MESH.edges_normal_sign, MESH.edges_non_unit_normals, MESH.edges_length] = normals(MESH);
MESH.gamma_ratio = gamma_ratio(MESH, PARAMETERS);
% The following value of tau depends on the tau-limiter. By default, we fix it to 1/gamma
tau = 1./MESH.gamma_ratio;
if (strcmp(PARAMETERS.DOMAIN.GEOMETRY, 'RECTANGLE') && ...
strcmp(PARAMETERS.MESH.GENERATION, 'NS2DDV') && ...
strcmp(PARAMETERS.FV.SCHEME, 'MUSCL') && ...
strcmp(PARAMETERS.FV.CELLS_DESIGN, 'SQUARES'))
[MESH.horizontal_edges, MESH.vertical_edges, MESH.volume, perimeter] = ...
mesh_squares_fv(MESH, PARAMETERS);
MESH.ciarlet_ratio = 2.;
MESH.coef_CFL = 1000; %%%% used only for testing renumbering techniques
%%%% used only for testing renumbering
for i=1:MESH.np1
MESH.coef_CFL = min(MESH.coef_CFL, ...
MESH.volume(i)/(perimeter(i)*(2.0+tau*MESH.gamma_ratio*MESH.ciarlet_ratio)));
end
MESH.volume(RENUMBERING.perm_p1) = MESH.volume;
MESH.volume_barycenter = MESH.p1;
else
[MESH.normal, MESH.volume, volume_P2, ...
perimeter, MESH.volume_barycenter, barycenter_patchEF] = mesh_stars_fv_part1(MESH);
% perimeter is local variable
% volume_P2 and barycenter_patchEF are only used for testing renumbering techniques
[MESH.boundary_edges, MESH.boundary_edges_cells, ...
boundary_edges_length, MESH.boundary_edges_normals] = mesh_stars_fv_part2(MESH);
% boundary_edges_length is only used for testing renumbering techniques
[upstream_cell, downstream_cell, MESH.adjacent_cells, MESH.ciarlet_ratio] ...
= mesh_stars_fv_part3(MESH);
% upstream_cell and downstream_cell are only used for testing renumbering techniques
MESH.coef_CFL = 1000; %%%% used only for testing renumbering techniques
for i=1:MESH.np1
%%%% used only for testing renumbering
MESH.coef_CFL = min(MESH.coef_CFL, ...
MESH.volume(i)/(perimeter(i)*(2.0+tau*MESH.gamma_ratio*MESH.ciarlet_ratio)));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Modify the descriptor adjacent_cells : this function locates %%%%%
%%%%% the upstream and downstream cells of each edge, %%%%%
%%%%% only for triangles which share a boundary node %%%%%
if strcmp(PARAMETERS.FV.FLUX_LIMITER,'TAULIM')
[MESH.adjacent_cells, MESH.flag_shadow_cells, MESH.corners] = ...
modify_adjacent_cells_border(MESH, PARAMETERS);
end
else
MESH.volume = mesh_stars_volume(MESH);
end
fprintf('done\n');
end
