function [fek,ind] = dexsets(deg,wam,leja,rect);

% extracts discrete extremal sets, namely approximate Fekete or Leja 
% interpolation points, from a 2d or 3d weakly-admissible mesh (wam)

% by Alvise Sommariva and Marco Vianello, University of Padova 
% October 30, 2014


% INPUT: 

% deg = interpolation degree

% wam = 2- or 3-column array of wam points coordinates, from which the 
% approximate Fekete or Leja points will be extracted

% leja=0: computes approximate Fekete points
% leja=1: computes approximate Leja points

% rect = 4- or 6-component vector such that the rectangle 
% [rect(1),rect(2)] x [rect(3),rect(4)] in 2d
% or [rect(1),rect(2)] x [rect(3),rect(4)] x [rect(5),rect(6)] in 3d
% contains the mesh
% if rect == [ ], it will be automatically tailored to the mesh 


% OUTPUT: 

% fek = 2- or 3-column array of approximate Fekete or Leja 
% points coordinates

% ind = indexes of the selected points in wam


% FUNCTION BODY

% setting rectangle if not defined 
if isempty(rect) 
rect(1)=min(wam(:,1));
rect(2)=max(wam(:,1));
rect(3)=min(wam(:,2));
rect(4)=max(wam(:,2));
if length(wam(1,:)) == 3
rect(5)=min(wam(:,3));
rect(6)=max(wam(:,3));
end;
end;

% dimension of the polynomial space
if length(wam(1,:)) == 2
dim=(deg+1)*(deg+2)/2;   
end;
if length(wam(1,:)) == 3
dim=(deg+1)*(deg+2)*(deg+3)/6;
end;

% approximate Fekete or Leja points extraction from wam  

% change of basis via 2 successive QR factorizations of the 
% Chebyshev-Vandermonde matrix 
[Q,R1,R2]=wamdop(deg,wam,rect);

% extracting the interpolation points 

if leja == 0
% approximate Fekete points and corresponding weights
% by QR factorization with column pivoting of Q'
w=Q'\ones(dim,1); 
ind=find(abs(w)>0);
fek=wam(ind,:);
else
% approximate Leja points and corresponding weights  
% by LU factorization with row pivoting of Q 
[L,U,perm]=lu(Q,'vector');
ind=perm(1:dim);
fek=wam(ind,:);
end;

end



