function [QuadTriangle] = QuadraturePointsTriangles(mesh)

% Quadrature nodes in the reference triangle:
[psi,eta,weights] = QuadLTri(4);

% Basis functions in the reference coordinate:
Basis = [1-psi(:)-eta(:), psi(:), eta(:)]; 

% Quadrature mesh.points in each triangle: 
% x_qp = x1 (1-\psi-\eat) + x2 \psi + x3 \eta
% y_qp = y1 (1-\psi-\eat) + y2 \psi + y3 \eta

QuadTriangle.xnodes = Basis * [ mesh.points(1,mesh.triangles(1,:)); ...
                                mesh.points(1,mesh.triangles(2,:)); ...
                                mesh.points(1,mesh.triangles(3,:)) ];
             
QuadTriangle.ynodes = Basis * [ mesh.points(2,mesh.triangles(1,:)); ...
                                mesh.points(2,mesh.triangles(2,:)); ...
                                mesh.points(2,mesh.triangles(3,:)) ];

QuadTriangle.weights = weights;

QuadTriangle.number = size(QuadTriangle.weights,2);

end
