%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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[st1b] = add_subtri_FE_P1B(MESH)

	% INPUT
	% MESH		An mesh with at least P1 and P1B nodes and triangles
	%
	% OUTPUTS
	% st1b		The descriptor of P1B subtriangles by P1B nodes

	% Define P1B subtriangles
	%%%%%%%%%%%%%%%%%%%%%%%%%
	st1b = zeros(3,3*MESH.nt1);
	for ke=1:MESH.nt1
		st1b(1,3*(ke-1)+1) = MESH.t1b(1,ke);
		st1b(2,3*(ke-1)+1) = MESH.t1b(2,ke);
		st1b(3,3*(ke-1)+1) = MESH.t1b(4,ke);

		st1b(1,3*(ke-1)+2) = MESH.t1b(2,ke);
		st1b(2,3*(ke-1)+2) = MESH.t1b(3,ke);
		st1b(3,3*(ke-1)+2) = MESH.t1b(4,ke);

		st1b(1,3*(ke-1)+3) = MESH.t1b(3,ke);
		st1b(2,3*(ke-1)+3) = MESH.t1b(1,ke);
		st1b(3,3*(ke-1)+3) = MESH.t1b(4,ke);
	end

end