read_2D.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[res] = read_2D(namefile, header)
% INPUT
% namefile The result file to be read
% header The diagnostic to be extracted
%
% OUTPUT
% res The extracted result
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 ~(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
% Reminder : in .mat files, the vector fields are saved as vectors. Their components are not separated
% In .h5 files, there are only scalar fields (the vector fields are built in .xdmf descriptors)
if usemat
A = load(namefile);
if exist(A.namemesh)
B = load(A.namemesh);
else
error(sprintf('Cannot find mesh file %s', A.namemesh));
end
else if usehdf5
% Extract the data
headers = h5info(namefile);
for i=1:numel(headers.Datasets)
name = headers.Datasets(i).Name;
value = h5read(namefile, strcat('/', name));
switch name
% Density (scalar)
case {'Density'}
A.rho = value;
A.NAME_RHO = 'Density';
case {'Density_exact'}
A.rhoex = value;
A.NAME_RHOEX = 'Density (exact)';
% Pressure (scalar)
case {'Pressure'}
A.p = value;
A.NAME_P = 'Pressure';
case {'Pressure_exact'}
A.pex = value;
A.NAME_PEX = 'Pressure (exact)';
% Components of pressure gradient
case {'Pressure_gradient_x'}
A.dxp = value;
case {'Pressure_gradient_x_exact'}
A.dxpex = value;
case {'Pressure_gradient_y'}
A.dyp = value;
case {'Pressure_gradient_y_exact'}
A.dypex = value;
% Components of velocity
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;
% Derivatives of the velocity (scalar)
case {'Velocity_x_dx'}
A.dxux = value;
A.NAME_DXUX = 'x-derivative of velocity (x-component)';
case {'Velocity_x_dx_exact'}
A.dxuexx = value;
A.NAME_DXUXEX = 'x-derivative of velocity (x-component) (exact)';
case {'Velocity_x_dy'}
A.dyux = value;
A.NAME_DYUX = 'y-derivative of velocity (x-component)';
case {'Velocity_x_dy_exact'}
A.dyuexx = value;
A.NAME_DYUXEX = 'y-derivative of velocity (x-component) (exact)';
case {'Velocity_y_dx'}
A.dxuy = value;
A.NAME_DXUY = 'x-derivative of velocity (y-component)';
case {'Velocity_y_dx_exact'}
A.dxuexy = value;
A.NAME_DXUYEX = 'x-derivative of velocity (y-component) (exact)';
case {'Velocity_y_dy'}
A.dyuy = value;
A.NAME_DYUY = 'y-derivative of velocity (y-component)';
case {'Velocity_y_dy_exact'}
A.dyuexy = value;
A.NAME_DYUYEX = 'y-derivative of velocity (y-component) (exact)';
% Shear rate (scalar)
case {'Shear_rate'}
A.shear = value;
A.NAME_SHEAR = 'Shear rate';
case {'Shear_rate_exact'}
A.shearex = value;
A.NAME_SHEAREX = 'Shear rate (exact)';
% Velocity magnitude (scalar)
case {'Velocity_magnitude'}
A.umod = value;
A.NAME_UMOD = 'Velocity magnitude';
case {'Velocity_magnitude_exact'}
A.umod = value;
A.NAME_UMODEX = 'Velocity magnitude (exact)';
% Vorticity (exact)
case {'Vorticity'}
A.omega = value;
A.NAME_OMEGA = 'Vorticity';
case {'Vorticity_exact'}
A.omegaex = value;
A.NAME_OMEGAEX = 'Vorticity (exact)';
% Stream function (scalar)
case {'Streamlines'}
A.psi = value;
A.NAME_PSI = 'Streamlines';
case {'Streamlines_exact'}
A.psiex = value;
A.NAME_PSIEX = 'Streamlines (exact)';
% Components of auxiliary velocity
case {'w_x'}
A.wx = value;
case {'w_y'}
A.wy = value;
% Auxiliary pressure (scalar)
case {'r'}
A.r = value;
A.NAME_R = 'r';
case {'Time'}
A.time = value;
end
end
if (isfield(A, 'ux') && isfield(A, 'uy'))
A.u = [A.ux', A.uy'];
A.NAME_U = 'Velocity';
end
if (isfield(A, 'uexx') && isfield(A, 'uexy'))
A.uex = [A.uexx', A.uexy'];
A.NAME_UEX = 'Velocity (exact)';
end
if (isfield(A, 'dxp') && isfield(A, 'dyp'))
A.gradp = [A.dxp', A.dyp'];
A.NAME_GRADP = 'Pressure gradient';
end
if (isfield(A, 'dxpex') && isfield(A, 'dypex'))
A.gradpex = [A.dxpex', A.dypex'];
A.NAME_GRADPEX = 'Pressure gradient (exact)';
end
if (isfield(A, 'wx') && isfield(A, 'wy'))
A.w = [A.wx', A.wy'];
A.NAME_W = 'w';
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
B.GEOMETRY = h5readatt(namemesh, '/', 'Geometry');
A.METHOD_U = h5readatt(namefile, '/', 'Method_velocity');
A.METHOD_P = h5readatt(namefile, '/', 'Method_pressure');
end
for i=1:numel(header)
if isfield(A, header{i})
res{i} = getfield(A, header{i});
elseif isfield(B.GRID, header{i})
res{i} = getfield(B.GRID, header{i});
elseif isfield(B, header{i})
res{i} = getfield(B, header{i});
else
error(sprintf('Cannot find field %s in input files %s and %s\n', header{i}, namefile, namemesh));
end
end
end
