plot_from_file.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[] = plot_from_file(namefile, XDISPLAY_LIST, LEVELS, XDISPLAY_METHOD, typeplot, nbrows, nbcols)
% INPUT
% namefile The input file to be used for generating the animation
% XDISPLAY_LIST List of diagnostics to be plotted
% LEVELS Isovalues associated to the diagnostics to be plotted
% XDISPLAY_METHOD Kind of visualization (single window or multiple window)
% typeplot Kind of scalar data representation (contour or pseudocolor)
% nbrows Number of rows in the visualization grid
% nbcols Number of columns in the visualization grid
close all;
cprec = 1e-8;
v = version('-release');
vnum = str2num(v(1:4));
usemat = strcmp(namefile(size(namefile,2)-3:size(namefile,2)), '.mat');
usehdf5 = strcmp(namefile(size(namefile,2)-2:size(namefile,2)), '.h5');
if ~(numel(XDISPLAY_LIST) == numel(LEVELS))
error("Problem in plot_from_file: Second and third arguments must be cell array with same shape");
end
if ~(usemat || usehdf5)
error('Only .mat or .h5 files can be handled');
elseif (usehdf5 && (vnum < 2011))
error('Cannot read .h5 file: please Matlab R2011a or newer version to proceed');
end
hasdensity = false;
if usemat
% Extract all data
A = load(namefile);
if isfield(A, 'rho')
hasdensity = true;
end
% Extract the mesh
if exist(A.namemesh)
B = load(A.namemesh);
else
error(sprintf('Cannot find mesh file %s', A.namemesh));
end
elseif usehdf5
% Extract all data
headers = h5info(namefile);
for i=1:numel(headers.Datasets)
name = headers.Datasets(i).Name;
value = h5read(namefile, strcat('/', name));
switch name
case {'Density'}
A.rho = value;
hasdensity = true;
case {'Density (exact)'}
A.rhoex = value;
case {'Pressure'}
A.p = value;
case {'Pressure_exact'}
A.pex = value;
case {'Velocity_x'}
A.ux = value;
case {'Velocity_x_exact'}
A.uexx = value;
case {'Velocity_y'}
A.uy = value;
case {'Velocity_y_exact'}
A.uexy = value;
case {'Velocity_x_dx'}
A.dxux = value;
case {'Velocity_x_dy'}
A.dyux = value;
case {'Velocity_y_dx'}
A.dxuy = value;
case {'Velocity_y_dy'}
A.dyuy = value;
case {'Velocity_x_dx_exact'}
A.dxuexx = value;
case {'Velocity_x_dy_exact'}
A.dyuexx = value;
case {'Velocity_y_dx_exact'}
A.dxuexy = value;
case {'Velocity_y_dy_exact'}
A.dyuexy = value;
case {'Shear_rate'}
A.shear = value;
case {'Shear_rate_exact'}
A.shearex = value;
case {'Velocity_magnitude'}
A.umod = value;
case {'Velocity_magnitude_exact'}
A.umod = value;
case {'Vorticity'}
A.omega = value;
case {'Vorticity_exact'}
A.omegaex = value;
case {'Streamlines'}
A.psi = value;
case {'Streamlines_exact'}
A.psiex = value;
case {'Pressure_gradient_x'}
A.dxp = value;
case {'Pressure_gradient_y'}
A.dyp = value;
case {'Pressure_gradient_x_exact'}
A.dxpex = value;
case {'Pressure_gradient_y_exact'}
A.dypex = value;
case {'Time'}
A.time = value;
end
end
% Complete with vector fields
if (isfield(A, 'ux') && isfield(A, 'uy'))
A.u = [A.ux', A.uy'];
end
if (isfield(A, 'uexx') && isfield(A, 'uexy'))
A.uex = [A.uexx', A.uexy'];
end
if (isfield(A, 'wx') && isfield(A, 'wy'))
A.w = [A.wx', A.wy'];
end
if (isfield(A, 'dxp') && isfield(A, 'dyp'))
A.gradp = [A.dxp', A.dyp'];
end
if (isfield(A, 'dxpex') && isfield(A, 'dypex'))
A.gradpex = [A.dxpex', A.dypex'];
end
% Extract the mesh
namemesh = h5readatt(namefile, '/', 'Mesh');
A.NAME_MESH = 'Mesh';
if exist(namemesh)
headers = h5info(namemesh);
else
error(sprintf('Cannot find mesh file %s', namemesh));
end
for i=1:numel(headers.Datasets)
name = headers.Datasets(i).Name;
value = h5read(namemesh, strcat('/', name));
switch name
case {'MESH_p1'}
B.GRID.p1 = value;
case {'MESH_t1'}
B.GRID.t1 = value;
B.GRID.t1 = B.GRID.t1 + ones(size(value));
case {'MESH_p1b'}
B.GRID.p1b = value;
case {'MESH_st1b'}
B.GRID.st1b = value;
B.GRID.st1b = B.GRID.st1b + ones(size(value));
case {'MESH_p2'}
B.GRID.p2 = value;
case {'MESH_st2'}
B.GRID.st2 = value;
B.GRID.st2 = B.GRID.st2 + ones(size(value));
end
end
% Extract the spatial discretizations
A.METHOD_U = h5readatt(namefile, '/', 'Method_velocity');
A.METHOD_P = h5readatt(namefile, '/', 'Method_pressure');
if hasdensity
A.METHOD_RHO = h5readatt(namefile, '/', 'Method_density');
end
end
TEST_ANALYTIC = isfield(A, 'pex');
XDISPLAY_CASES = '''MESH'', ''SUBTRI_U'', ';
XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''P'', ''P_DX'', ''P_DY'', ''GRAD_P'', ');
XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''U_VECTOR'', ''UX'', ''UY'', ');
XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''VORTIC'', ''STREAM'', ''U_MODULE'', ');
XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''UX_DX'', ''UX_DY'', ''UY_DX'', ''UY_DY'', ''SHEAR'', ');
XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''W_VECTOR'', ''WX'', ''WY'', ''R''');
if hasdensity
XDISPLAY_CASES = strcat(XDISPLAY_CASES, ', ''RHO''');
end
if TEST_ANALYTIC
XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''P_EX'', ''P_DX_EX'', ''P_DY_EX'', ''GRAD_P_EX'', ');
XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''U_VECTOR_EX'', ''UX_EX'', ''UY_EX'', ');
XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''VORTIC_EX'', ''STREAM_EX'', ''U_MODULE_EX'', ');
XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''UX_DX_EX'', ''UX_DY_EX'', ''UY_DX_EX'', ''UY_DY_EX'', ');
XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''SHEAR_EX''');
if hasdensity
XDISPLAY_CASES = strcat(XDISPLAY_CASES, ', ''RHO_EX''');
end
end
% Prepare the isovalues for scalar diagnostics
ind_rho = [];
ind_p = [];
ind_dxp = [];
ind_dyp = [];
ind_ux = [];
ind_uy = [];
ind_omega = [];
ind_psi = [];
ind_shear = [];
ind_umod = [];
ind_dxux = [];
ind_dyux = [];
ind_dxuy = [];
ind_dyuy = [];
for i=1:numel(XDISPLAY_LIST)
switch XDISPLAY_LIST{i}
case {'MESH', 'SUBTRI_U', 'U_VECTOR', 'U_VECTOR_EX', 'W_VECTOR', 'GRAD_P'}
% Nothing to do
case {'RHO', 'RHO_EX'}
ind_rho = [ind_rho, i];
if TEST_ANALYTIC
levels_rho = compute_levels_from_file(A.rho, A.rhoex, ind_rho, LEVELS);
else
levels_rho = compute_levels_from_file(A.rho, [], ind_rho, LEVELS);
end
case {'P', 'P_EX'}
ind_p = [ind_p, i];
if TEST_ANALYTIC
levels_p = compute_levels_from_file(A.p, A.pex, ind_p, LEVELS);
else
levels_p = compute_levels_from_file(A.p, [], ind_p, LEVELS);
end
case {'P_DX', 'P_DX_EX'}
ind_dxp = [ind_dxp, i];
if TEST_ANALYTIC
levels_dxp = compute_levels_from_file(A.dxp, A.dxpex, ind_dxp, LEVELS);
else
levels_dxp = compute_levels_from_file(A.dxp, [], ind_dxp, LEVELS);
end
case {'P_DY', 'P_DY_EX'}
ind_dyp = [ind_dyp, i];
if TEST_ANALYTIC
levels_dyp = compute_levels_from_file(A.dyp, A.dypex, ind_dyp, LEVELS);
else
levels_dyp = compute_levels_from_file(A.dyp, [], ind_dyp, LEVELS);
end
case {'UX', 'UX_EX'}
ind_ux = [ind_ux, i];
if TEST_ANALYTIC
levels_ux = compute_levels_from_file(A.ux, A.uxex, ind_ux, LEVELS);
else
levels_ux = compute_levels_from_file(A.ux, [], ind_ux, LEVELS);
end
case {'UY', 'UY_EX'}
ind_uy = [ind_uy, i];
if TEST_ANALYTIC
levels_uy = compute_levels_from_file(A.uy, A.uyex, ind_uy, LEVELS);
else
levels_uy = compute_levels_from_file(A.uy, [], ind_uy, LEVELS);
end
case {'VORTIC', 'VORTIC_EX'}
ind_omega = [ind_omega, i];
if TEST_ANALYTIC
levels_omega = compute_levels_from_file(A.omega, A.omegaex, ind_omega, LEVELS);
else
levels_omega = compute_levels_from_file(A.omega, [], ind_omega, LEVELS);
end
case {'STREAM', 'STREAM_EX'}
ind_psi = [ind_psi, i];
if TEST_ANALYTIC
levels_psi = compute_levels_from_file(A.psi, A.psiex, ind_psi, LEVELS);
else
levels_psi = compute_levels_from_file(A.psi, [], ind_psi, LEVELS);
end
case {'U_MODULE', 'U_MODULE_EX'}
ind_umod = [ind_umod, i];
if TEST_ANALYTIC
levels_umod = compute_levels_from_file(A.umod, A.umodex, ind_umod, LEVELS);
else
levels_umod = compute_levels_from_file(A.umod, [], ind_umod, LEVELS);
end
case {'UX_DX', 'UX_DX_EX'}
ind_dxux = [ind_dxux, i];
if TEST_ANALYTIC
levels_dxux = compute_levels_from_file(A.dxux, A.dxuxex, ind_dxux, LEVELS);
else
levels_dxux = compute_levels_from_file(A.dxux, [], ind_dxux, LEVELS);
end
case {'UX_DY', 'UX_DY_EX'}
ind_dyux = [ind_dyux, i];
if TEST_ANALYTIC
levels_dyux = compute_levels_from_file(A.dyux, A.dyuxex, ind_dyux, LEVELS);
else
levels_dyux = compute_levels_from_file(A.dyux, [], ind_dyux, LEVELS);
end
case {'UY_DX', 'UY_DX_EX'}
ind_dxuy = [ind_dxuy, i];
if TEST_ANALYTIC
levels_dxuy = compute_levels_from_file(A.dxuy, A.dxuyex, ind_dxuy, LEVELS);
else
levels_dxuy = compute_levels_from_file(A.dxuy, [], ind_dxuy, LEVELS);
end
case {'UY_DY', 'UY_DY_EX'}
ind_dyuy = [ind_dyuy, i];
if TEST_ANALYTIC
levels_dyuy = compute_levels_from_file(A.dyuy, A.dyuyex, ind_dyuy, LEVELS);
else
levels_dyuy = compute_levels_from_file(A.dyuy, [], ind_dyuy, LEVELS);
end
case {'SHEAR', 'SHEAR_EX'}
ind_shear = [ind_shear, i];
if TEST_ANALYTIC
levels_shear = compute_levels_from_file(A.shear, A.shearex, ind_shear, LEVELS);
else
levels_shear = compute_levels_from_file(A.shear, [], ind_shear, LEVELS);
end
case {'WX'}
levels_wx = compute_levels_from_file(A.wx, [], [i], LEVELS);
case {'WY'}
levels_wy = compute_levels_from_file(A.wy, [], [i], LEVELS);
case {'R'}
levels_r = compute_levels_from_file(A.r, [], [i], LEVELS);
otherwise
error(sprintf('Cannot plot the diagnostic %s\nOnly the following diagnostics are allowed: %s', ...
XDISPLAY_LIST{i}, XDISPLAY_CASES));
end
end
switch XDISPLAY_METHOD
case {'SINGLE_FIGURE'}
% Start the visualization interface
if ((~exist('nbcols')) || (~exist('nbrows')))
error(sprintf('If SINGLE_FIGURE method is selected, you must provide a number of rows and a number of columns\n'));
end
HANDLES = start_figs(XDISPLAY_LIST, XDISPLAY_METHOD, nbrows, nbcols);
% Add the time
figure(HANDLES.h);
globtitle = sprintf('Time t = %e', A.time);
meshtitle = HANDLES.meshtitle;
subtriutitle = sprintf('%s %s', A.METHOD_U, HANDLES.subtriutitle);
uvecttitle = HANDLES.uvecttitle;
uvectextitle = HANDLES.uvectextitle;
uxtitle = HANDLES.uxtitle;
uxextitle = HANDLES.uxextitle;
uytitle = HANDLES.uytitle;
uyextitle = HANDLES.uyextitle;
dxuxtitle = HANDLES.dxuxtitle;
dxuxextitle = HANDLES.dxuxextitle;
dyuxtitle = HANDLES.dyuxtitle;
dyuxextitle = HANDLES.dyuxextitle;
dxuytitle = HANDLES.dxuytitle;
dxuyextitle = HANDLES.dxuyextitle;
dyuytitle = HANDLES.dyuytitle;
dyuyextitle = HANDLES.dyuyextitle;
sheartitle = HANDLES.sheartitle;
shearextitle = HANDLES.shearextitle;
ptitle = HANDLES.ptitle;
pextitle = HANDLES.pextitle;
dxptitle = HANDLES.dxptitle;
dxpextitle = HANDLES.dxpextitle;
dyptitle = HANDLES.dyptitle;
dypextitle = HANDLES.dypextitle;
gradptitle = HANDLES.gradptitle;
gradpextitle = HANDLES.gradpextitle;
omegatitle = HANDLES.omegatitle;
omegaextitle = HANDLES.omegaextitle;
psititle = HANDLES.psititle;
psiextitle = HANDLES.psiextitle;
umodtitle = HANDLES.umodtitle;
umodextitle = HANDLES.umodextitle;
rhotitle = HANDLES.rhotitle;
rhoextitle = HANDLES.rhoextitle;
rtitle = HANDLES.rtitle;
wvecttitle = HANDLES.wvecttitle;
wxtitle = HANDLES.wxtitle;
wytitle = HANDLES.wytitle;
case {'MULTIPLE_FIGURE'}
% Start the visualization interface
HANDLES = start_figs(XDISPLAY_LIST, XDISPLAY_METHOD);
meshtitle = sprintf('%s at time = %e', HANDLES.meshtitle, A.time);
subtriutitle = sprintf('%s %s at time = %e', A.METHOD_U, HANDLES.subtriutitle, A.time);
uvecttitle = sprintf('%s at time = %e', HANDLES.uvecttitle, A.time);
uvectextitle = sprintf('%s at time = %e', HANDLES.uvectextitle, A.time);
uxtitle = sprintf('%s at time = %e', HANDLES.uxtitle, A.time);
uxextitle = sprintf('%s at time = %e', HANDLES.uxextitle, A.time);
uytitle = sprintf('%s at time = %e', HANDLES.uytitle, A.time);
uyextitle = sprintf('%s at time = %e', HANDLES.uyextitle, A.time);
dxuxtitle = sprintf('%s at time = %e', HANDLES.dxuxtitle, A.time);
dxuxextitle = sprintf('%s at time = %e', HANDLES.dxuxextitle, A.time);
dyuxtitle = sprintf('%s at time = %e', HANDLES.dyuxtitle, A.time);
dyuxextitle = sprintf('%s at time = %e', HANDLES.dyuxextitle, A.time);
dxuytitle = sprintf('%s at time = %e', HANDLES.dxuytitle, A.time);
dxuyextitle = sprintf('%s at time = %e', HANDLES.dxuyextitle, A.time);
dyuytitle = sprintf('%s at time = %e', HANDLES.dyuytitle, A.time);
dyuyextitle = sprintf('%s at time = %e', HANDLES.dyuyextitle, A.time);
sheartitle = sprintf('%s at time = %e', HANDLES.sheartitle, A.time);
shearextitle = sprintf('%s at time = %e', HANDLES.shearextitle, A.time);
ptitle = sprintf('%s at time = %e', HANDLES.ptitle, A.time);
pextitle = sprintf('%s at time = %e', HANDLES.pextitle, A.time);
dxptitle = sprintf('%s at time = %e', HANDLES.dxptitle, A.time);
dxpextitle = sprintf('%s at time = %e', HANDLES.dxpextitle, A.time);
dyptitle = sprintf('%s at time = %e', HANDLES.dyptitle, A.time);
dypextitle = sprintf('%s at time = %e', HANDLES.dypextitle, A.time);
gradptitle = sprintf('%s at time = %e', HANDLES.gradptitle, A.time);
gradpextitle = sprintf('%s at time = %e', HANDLES.gradpextitle, A.time);
omegatitle = sprintf('%s at time = %e', HANDLES.omegatitle, A.time);
omegaextitle = sprintf('%s at time = %e', HANDLES.omegaextitle, A.time);
psititle = sprintf('%s at time = %e', HANDLES.psititle, A.time);
psiextitle = sprintf('%s at time = %e', HANDLES.psiextitle, A.time);
umodtitle = sprintf('%s at time = %e', HANDLES.umodtitle, A.time);
umodextitle = sprintf('%s at time = %e', HANDLES.umodextitle, A.time);
rhotitle = sprintf('%s at time = %e', HANDLES.rhotitle, A.time);
rhoextitle = sprintf('%s at time = %e', HANDLES.rhoextitle, A.time);
rtitle = sprintf('%s at time = %e', HANDLES.rtitle, A.time);
wvecttitle = sprintf('%s at time = %e', HANDLES.wvecttitle, A.time);
wxtitle = sprintf('%s at time = %e', HANDLES.wxtitle, A.time);
wytitle = sprintf('%s at time = %e', HANDLES.wytitle, A.time);
otherwise
error(sprintf('Problem while using plot_from_file: wrong value of PARAMETERS.OUTPUT.XDISPLAY\n'));
end
for i=1:numel(XDISPLAY_LIST)
% Identify handle, levels, title, data
switch XDISPLAY_LIST{i}
case {'MESH'}
h = HANDLES.hmesh;
ti = meshtitle;
case {'SUBTRI_U'}
h = HANDLES.hsubtriu;
ti = subtriutitle;
case {'RHO'}
h = HANDLES.hrho;
levels = levels_rho;
ti = rhotitle;
data = A.rho;
case {'RHO_EX'}
h = HANDLES.hrhoex;
levels = levels_rho;
ti = rhoextitle;
data = A.rhoex;
case {'P'}
h = HANDLES.hp;
levels = levels_p;
ti = ptitle;
data = A.p;
case {'P_EX'}
h = HANDLES.hpex;
levels = levels_p;
ti = pextitle;
data = A.pex;
case {'P_DX'}
h = HANDLES.hdxp;
levels = levels_dxp;
ti = dxptitle;
data = A.dxp;
case {'P_DX_EX'}
h = HANDLES.hdxpex;
levels = levels_dxpex;
ti = dxpextitle;
data = A.dxpex;
case {'P_DY'}
h = HANDLES.hdyp;
levels = levels_dyp;
ti = dyptitle;
data = A.dyp;
case {'P_DY_EX'}
h = HANDLES.hdypex;
levels = levels_dypex;
ti = dypextitle;
data = A.dypex;
case {'GRAD_P'}
h = HANDLES.hgradp;
ti = gradptitle;
data = A.gradp;
case {'GRAD_P_EX'}
h = HANDLES.hgradpex;
ti = gradpextitle;
data = A.gradpex;
case {'U_VECTOR'}
h = HANDLES.huvect;
ti = uvecttitle;
data = A.u;
case {'U_VECTOR_EX'}
h = HANDLES.huvectex;
ti = uvectextitle;
data = A.uex;
case {'UX'}
h = HANDLES.hux;
levels = levels_ux;
ti = uxtitle;
data = A.ux;
case {'UX_EX'}
h = HANDLES.huxex;
levels = levels_ux;
ti = uxextitle;
data = A.uxex;
case {'UY'}
h = HANDLES.huy;
levels = levels_uy;
ti = uytitle;
data = A.uy;
case {'UY_EX'}
h = HANDLES.huy_ex;
levels = levels_uy;
ti = uyextitle;
data = A.uyex;
case {'UX_DX'}
h = HANDLES.hdxux;
levels = levels_dxux;
ti = dxuxtitle;
data = A.dxux;
case {'UX_DX_EX'}
h = HANDLES.hdxuxex;
levels = levels_dxux;
ti = dxuxextitle;
data = A.dxuxex;
case {'UX_DY'}
h = HANDLES.hdyux;
levels = levels_dyux;
ti = dyuxtitle;
data = A.dyux;
case {'UX_DY_EX'}
h = HANDLES.hdyuxex;
levels = levels_dyux;
ti = dyuxextitle;
data = A.dyuxex;
case {'UY_DX'}
h = HANDLES.hdxuy;
levels = levels_dxuy;
ti = dxuytitle;
data = A.dxuy;
case {'UY_DX_EX'}
h = HANDLES.hdxuyex;
levels = levels_dxuy;
ti = dxuyextitle;
data = A.dxuyex;
case {'UY_DY'}
h = HANDLES.hdyuy;
levels = levels_dyuy;
ti = dyuytitle;
data = A.dyuy;
case {'UY_DY_EX'}
h = HANDLES.hdyuyex;
levels = levels_dyuy;
ti = dyuyextitle;
data = A.dyuyex;
case {'SHEAR'}
h = HANDLES.hshear;
levels = levels_shear;
ti = sheartitle;
data = A.shear;
case {'SHEAR_EX'}
h = HANDLES.hshearex;
levels = levels_shear;
ti = shearextitle;
data = A.shearex;
case {'VORTIC'}
h = HANDLES.homega;
levels = levels_omega;
ti = omegatitle;
data = A.omega;
case {'VORTIC_EX'}
h = HANDLES.homegaex;
levels = levels_omega;
ti = omegaextitle;
data = A.omegaex;
case {'STREAM'}
h = HANDLES.hpsi;
levels = levels_psi;
ti = psititle;
data = A.psi;
case {'STREAM_EX'}
h = HANDLES.hpsiex;
levels = levels_psi;
ti = psiextitle;
data = A.psiex;
case {'U_MODULE'}
h = HANDLES.humod;
levels = levels_umod;
ti = umodtitle;
data = A.umod;
case {'U_MODULE_EX'}
h = HANDLES.humodex;
levels = levels_umod;
ti = umodextitle;
data = A.umodex;
case {'W_VECTOR'}
h = HANDLES.hwvect;
ti = wvecttitle;
data = A.w;
case {'WX'}
h = HANDLES.hwx;
levels = levels_wx;
ti = wxtitle;
data = A.wx;
case {'WY'}
h = HANDLES.hwy;
levels = levels_uy;
ti = wytitle;
data = A.wy;
case {'R'}
h = HANDLES.hr;
levels = levels_r;
ti = rtitle;
data = A.r;
otherwise
error(sprintf('Cannot plot the diagnostic %s\nOnly the following diagnostics are allowed: %s', ...
XDISPLAY_LIST{i}, XDISPLAY_CASES));
end
% Identify the mesh
switch XDISPLAY_LIST{i}
case {'MESH', 'R', 'P', 'P_EX', 'P_DX', 'P_DX_EX', 'P_DY', 'P_DY_EX', ...
'VORTIC', 'VORTIC_EX', 'STREAM', 'STREAM_EX', ...
'UX_DX', 'UX_DX_EX', 'UX_DY', 'UX_DY_EX', 'UY_DX', 'UY_DX_EX', ...
'UY_DY', 'UY_DY_EX', 'SHEAR', 'SHEAR_EX', 'RHO', 'RHO_EX'}
p = B.GRID.p1;
t = B.GRID.t1(1:3,:);
case {'GRAD_P', 'GRAD_P_EX'}
p = zeros(2,size(B.GRID.t1,2));
for ke=1:size(p,2)
p(:,ke) = (B.GRID.p1(:,B.GRID.t1(1,ke)) + B.GRID.p1(:,B.GRID.t1(2,ke)) + ...
B.GRID.p1(:,B.GRID.t1(3,ke)))/3.;
end
case {'SUBTRI_U', 'U_VECTOR', 'U_VECTOR_EX', 'UX', 'UX_EX', 'UY', 'UY_EX', ...
'U_MODULE', 'U_MODULE_EX', 'W_VECTOR', 'WX', 'WY'}
switch A.METHOD_U
case {'P1B'}
p = B.GRID.p1b;
t = B.GRID.st1b;
case {'P2'}
p = B.GRID.p2;
t = B.GRID.st2;
otherwise
error('Unknown FE discretization');
end
otherwise
error(sprintf('Cannot plot the diagnostic %s\nOnly the following diagnostics are allowed: %s', ...
XDISPLAY_LIST{i}, XDISPLAY_CASES));
end
% Plot data
switch_plot(h, XDISPLAY_METHOD);
switch XDISPLAY_LIST{i}
case {'MESH', 'SUBTRI_U'}
% Plot a mesh
plot_mesh(p, t, ti, h);
case {'R', 'RHO', 'RHO_EX', 'P', 'P_EX', 'VORTIC', 'VORTIC_EX', 'STREAM', 'STREAM_EX', ...
'UX_DX', 'UX_DY', 'UY_DX', 'UY_DY', 'SHEAR', ...
'UX_DX_EX', 'UX_DY_EX', 'UY_DX_EX', 'UY_DY_EX', 'SHEAR_EX', ...
'UX', 'UY', 'U_MODULE', 'UX_EX', 'UY_EX', 'U_MODULE_EX', 'WX', 'WY'}
% Plot a node-centered scalar data
switch typeplot
case {'contour'}
plot_contour(data, p, t, levels, ti, h);
case {'pseudocolor'}
plot_color(data, p, t, ti, h);
otherwise
error('Only ''contour'' or ''pseudocolor'' plot type is allowed');
end
case {'P_DX', 'P_DX_EX', 'P_DY', 'P_DY_EX'}
% Plot a cell-centered scalar data
plot_color(data, p, t, ti, h);
case {'U_VECTOR', 'U_VECTOR_EX', 'W_VECTOR', 'GRAD_P', 'GRAD_P_EX'}
% Plot a vector field
plot_vector(data, p, ti, h);
otherwise
error(sprintf('Cannot plot the diagnostic %s\nOnly the following diagnostics are allowed: %s', ...
XDISPLAY_LIST{i}, XDISPLAY_CASES));
end
if strcmp(XDISPLAY_METHOD, 'SINGLE_FIGURE')
suplabel(globtitle, 't');
drawnow;
end
end
%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary functions %
%%%%%%%%%%%%%%%%%%%%%%%
function[newlevels] = compute_levels_from_file(data, dataex, ind, LEVELS)
cprec = 1e-8;
newlevels = [];
if (numel(dataex) > 0)
% We wait for the reading of both data and exact data
if (numel(ind) == 2)
if (numel(LEVELS{ind(1)}) == numel(LEVELS{ind(2)}))
if (max(abs(LEVELS{ind(1)}-LEVELS{ind(2)}))/max(abs(LEVELS{ind(1)})) < cprec)
newlevels = LEVELS{ind(1)};
else
% Need to correct the isovalues
isomin = min(min([data,dataex]));
isomax = max(max([data,dataex]));
if (numel(LEVELS{ind(1)}) > 0)
newlevels = linspace(isomin,isomax,numel(LEVELS{ind(1)}));
else
newlevels = linspace(isomin,isomax,7);
end
end
else
% Need to correct the isovalues
isomin = min(min([data,dataex]));
isomax = max(max([data,dataex]));
newlevels = linspace(isomin,isomax,7);
end
end
else
if (numel(LEVELS{ind(1)}) > 0)
newlevels = LEVELS{ind(1)};
else
% Need to compute the isovalues
isomin = min(min(data));
isomax = max(max(data));
newlevels = linspace(isomin,isomax,7);
end
end
end
end
