%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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[] = save_mesh(time, MESH, namefile, PARAMETERS)

	% INPUT
	% time			Time value
	% MESH			Mesh metadata
	% namefile		Name of output mesh file
	% PARAMETERS	List of simulation parameters

	v = version('-release');
	vnum = str2num(v(1:4));
	usemat = (vnum < 2011);
	usehdf5 = ~usemat;

	METHOD_U = PARAMETERS.FE.TYPE;
   	METHOD_P = 'P1';
	Geometry = PARAMETERS.DOMAIN.GEOMETRY;

	if usemat
	    save(namefile, 'time');
    	NAME_MESH = 'Mesh';

    	save(namefile, 'NAME_MESH', 'METHOD_U', 'METHOD_P', 'Geometry', '-append');
    	GRID.p1 = MESH.p1;
    	GRID.t1 = MESH.t1;
    	if strcmp(METHOD_U, 'P1B')
    	    GRID.p1b = MESH.p1b;
    	    GRID.st1b = MESH.st1b;
    	end
    	if strcmp(METHOD_U, 'P2')
    	    GRID.p2 = MESH.p2;
    	    GRID.st2 = MESH.st2;
    	end
    	save(namefile, 'GRID', '-append');
	end

	if usehdf5
		nameh5 = strcat(namefile(1:(size(namefile,2)-3)), 'h5');
		unix(sprintf('rm %s 2>/dev/null', nameh5));

		h5create(nameh5, '/MESH_p1', size(MESH.p1), 'Datatype', 'double');
		h5write(nameh5, '/MESH_p1', MESH.p1);
		h5create(nameh5, '/MESH_t1', [3,MESH.nt1], 'Datatype', 'double');
		h5write(nameh5, '/MESH_t1', (MESH.t1(1:3,1:MESH.nt1)-ones(3,MESH.nt1)));

		 if strcmp(METHOD_U, 'P1B')
			h5create(nameh5, '/MESH_p1b', size(MESH.p1b), 'Datatype', 'double');
			h5write(nameh5, '/MESH_p1b', MESH.p1b);
			h5create(nameh5, '/MESH_st1b', size(MESH.st1b), 'Datatype', 'double');
			h5write(nameh5, '/MESH_st1b', MESH.st1b-ones(size(MESH.st1b)));
    	end
    	if strcmp(METHOD_U, 'P2')
			h5create(nameh5, '/MESH_p2', size(MESH.p2), 'Datatype', 'double');
			h5write(nameh5, '/MESH_p2', MESH.p2);
			h5create(nameh5, '/MESH_st2', size(MESH.st2), 'Datatype', 'double');
			h5write(nameh5, '/MESH_st2', MESH.st2-ones(size(MESH.st2)));
    	end

		h5writeatt(nameh5, '/', 'Geometry', Geometry);
	end
end