%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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[p2, np2, t2, nt2] = add_mdmesh_FE_P2(MESH)

	% INPUT
	% MESH		An mesh with at least P1 nodes and triangles
	%
	% OUTPUTS
	% p2		The P2 nodes
	% np2		The number of P2 nodes
	% t2		The triangle descriptor by P2 nodes
	% nt2		The number of P2 triangles (same as the number of P1 triangles)

	% Define the P2 list of nodes
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Define the center of each edge of the mesh
	pedges = zeros(2,MESH.ne1full);
	pedges(1,:) = (MESH.p1(1,MESH.e1full(1,:)) + MESH.p1(1,MESH.e1full(2,:)))/2.0; 
	pedges(2,:) = (MESH.p1(2,MESH.e1full(1,:)) + MESH.p1(2,MESH.e1full(2,:)))/2.0;
	% Add it to the MESH structure
	p2 = [MESH.p1 pedges];
	np2 = MESH.np1 + MESH.ne1full;

	% Define the P2 list of triangles
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Prepare the indices of the edge middle points
	tmiddle = abs(MESH.t1e) + MESH.np1*ones(3,MESH.nt1);
	% Enrich the P1 triangle descriptor with the the indices of edge centers sorted by triangles as P2 nodes
	t2 = [MESH.t1(1:3,:); tmiddle];
	nt2 = MESH.nt1;

	clear pedges;
	clear tmiddle;
end