
function [nodes,weights]=define_rule_simplex(ade)

% Cubature rule on the unit simplex 
% * with vertices (0,0), (1,0), (1,1),
% * with ADE equal to ade.

% Gaussian-Jacobi rule.
m=ceil((ade+1)/2);
ab_GJ=r_jacobi(m,0,1); xw_GJ=gauss(m,ab_GJ); % Gauss-Jacobi
x_GJ=(xw_GJ(:,1)+1)/2; w_GJ=xw_GJ(:,2);

% Gaussian-Legendre rule.
ab_GL=r_jacobi(m,0,0); xw_GL=gauss(m,ab_GL); % Gauss-Legendre
x_GL=(xw_GL(:,1)+1)/2; w_GL=xw_GL(:,2);

% Define tensorial rule
[x_mat_GJ,x_mat_GL]=meshgrid(x_GJ,x_GL);
X_mat=x_mat_GJ; Y_mat=x_mat_GJ.*x_mat_GL;

[w_mat_GJ,w_mat_GL]=meshgrid(w_GJ,w_GL);
W_mat=(1/8)*w_mat_GJ.*w_mat_GL;

nodes=[X_mat(:) Y_mat(:)]; 
weights=W_mat(:);