%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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[p1b, np1b, t1b, nt1b] = add_mdmesh_FE_P1B(MESH)

	% INPUT
	% MESH		An mesh with at least P1 nodes and triangles
	%
	% OUTPUTS
	% p1b		The P1B nodes
	% np1b		The number of P1B nodes
	% t1b		The triangle descriptor by P1B nodes
	% nt1b		The number of P1B triangles (same as the number of P1 triangles)

	% Define the P1B list of nodes
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Define the gravity center of each P1 triangle
	pbary = zeros(2,MESH.nt1);
	pbary(1,:) = (MESH.p1(1,MESH.t1(1,:)) + MESH.p1(1,MESH.t1(2,:)) + MESH.p1(1,MESH.t1(3,:)))/3.;
	pbary(2,:) = (MESH.p1(2,MESH.t1(1,:)) + MESH.p1(2,MESH.t1(2,:)) + MESH.p1(2,MESH.t1(3,:)))/3.;
	% Add it to the MESH structure
	p1b = [MESH.p1 pbary];
	np1b = MESH.np1 + MESH.nt1;

	% Define the P1B list of triangles
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Prepare the indices of the triangle centers in the P1B mesh
	tbary = (1:MESH.nt1) + MESH.np1*ones(1,MESH.nt1);
	% Enrich the P1 triangle descriptor with the the indices of triangle centers as P1B nodes
	t1b = [MESH.t1(1:3,:) ; tbary];
	nt1b = MESH.nt1;

	clear pbary;
	clear tbary;
end