Raw File
movie_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[] = movie_from_file(inputs, XDISPLAY_LIST, LEVELS, XDISPLAY_METHOD, typeplot, nbrows, nbcols)

	% INPUT
	% inputs				List of input files 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));

	errmsg = sprintf('Problem: first argument of movie_from_file must be a cell array structured as follows:\n');
	errmsg = strcat(errmsg, '{''/PATH/TO/RESULTS/prefix_*.h5''}\n');
	errmsg = strcat(errmsg, 'or\n');
	errmsg = strcat(errmsg, '{''/PATH/TO/RESULTS/prefix_*.h5'', index_start, index_end}\n');
	errmsg = strcat(errmsg, 'or\n');
	errmsg = strcat(errmsg, '{''/PATH/TO/RESULTS/prefix_*.h5'', index_start, index_end, step}\n');

	switch numel(inputs)
		case {1}
			pattern = inputs{1};
			step = 1;

			dirlistfile = dir(pattern);
			if (numel(dirlistfile) == 0)
				error(sprintf('Cannot find file with pattern %s\n', pattern));
			end

			fullpath = dirlistfile(1).folder;
			splitpattern = strsplit(pattern, '/');
			patternfile = splitpattern{numel(splitpattern)};
			splitpatternfile = strsplit(patternfile, '*');
			prefix = splitpatternfile{1};
			ext = splitpatternfile{2};

			indarray = [];
			for i=1:numel(dirlistfile)
				splitname = strsplit(dirlistfile(i).name, prefix);
				suffix = splitname{2};
				splitsuffix = strsplit(suffix, ext);
				indarray = [indarray, str2num(splitsuffix{1})];
			end
			indarray = sort(indarray);

			listfile = cell(1,numel(indarray));
			for i=1:numel(indarray)
				listfile{i} = strcat(fullpath, '/', prefix, num2str(indarray(i)), ext);
				lstest = (numel(dir(listfile{i})) > 0);
				if ~lstest
					error(sprintf('Problem with movie_from_file: cannot find file %s\n', listfile{i}));
				end
			end

		case {3}
			pattern = inputs{1};
			istart = inputs{2};
			iend = inputs{3};
			step = 1;

			dirlistfile = dir(pattern);
			if (numel(dirlistfile) == 0)
				error(sprintf('Cannot find file with pattern %s\n', pattern));
			end

			fullpath = dirlistfile(1).folder;
			splitpattern = strsplit(pattern, '/');
			patternfile = splitpattern{numel(splitpattern)};
			splitpatternfile = strsplit(patternfile, '*');
			prefix = splitpatternfile{1};
			ext = splitpatternfile{2};
			
			indarray = istart:iend;

			listfile = cell(1,numel(indarray));
			for i=1:numel(indarray)
				listfile{i} = strcat(fullpath, '/', prefix, num2str(indarray(i)), ext);
				lstest = (numel(dir(listfile{i})) > 0);
				if ~lstest
					error(sprintf('Problem with movie_from_file: cannot find file %s\n', listfile{i}));
				end
			end

		case {4}
			pattern = inputs{1};
			istart = inputs{2};
			iend = inputs{3};
			step = inputs{4};

			dirlistfile = dir(pattern);
			if (numel(dirlistfile) == 0)
				error(sprintf('Cannot find file with pattern %s\n', pattern));
			end

			fullpath = dirlistfile(1).folder;
			splitpattern = strsplit(pattern, '/');
			patternfile = splitpattern{numel(splitpattern)};
			splitpatternfile = strsplit(patternfile, '*');
			prefix = splitpatternfile{1};
			ext = splitpatternfile{2};

			indarray = istart:step:iend;

			listfile = cell(1,numel(indarray));
			for i=1:numel(indarray)
				listfile{i} = strcat(fullpath, '/', prefix, num2str(indarray(i)), ext);
				lstest = (numel(dir(listfile{i})) > 0);
				if ~lstest
					error(sprintf('Problem with movie_from_file: cannot find file %s\n', listfile{i}));
				end
			end

		otherwise
			error(errmsg);
	end

	usemat = strcmp(ext, '.mat');
	usehdf5 = strcmp(ext, '.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;

	for ifile=1:numel(listfile)
		namefile = listfile{ifile};

		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
				if (ifile == 1)
					HANDLES = start_figs(XDISPLAY_LIST, XDISPLAY_METHOD, nbrows, nbcols);
				end
				% 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
				if (ifile == 1)
					HANDLES = start_figs(XDISPLAY_LIST, XDISPLAY_METHOD);
				end

				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

		% Pauses the visualization during 0.5 seconds
		pause(0.1)
	end % End of list of input files


	%%%%%%%%%%%%%%%%%%%%%%%
	% 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
back to top