function [XIout,XIindexout,nsp] = FEconstruct(xknots,fval,order,ngauss,di,dj,constr,varargin)
%function [XI,XIindex,nsp] = FEconstruct(xknots,fval,sporder,ngauss,di,dj,constr,xknots_core)
%
% Each element of a FE matrix has the form int(Li*Lj*fval*p(t)) dx
% p(t) is the only time-varying part.
%
% This function constructs matrix XI such that
% the elements on the diagonals of the FE matrix are Pdiag=XI*p(xgauss);
% where p is time-varying part of the integrand evaluated on the gauss
% points xgauss.
% XIindex contains the linear indices into FE matrix P, pointing to the
% location where each value of Pdiag should go, i.e. P(XIindex) = XI*p;
% nsp is the number of splines
%
%  If called with dj=[]; then returns matrix for the source term
%    S = XI*s(xgauss);
%  if fval=scalar then the fixed part of the integrand is assumed not to
%    depend on rho.
%
%  INPUTS
%   xknots: spline knots
%   fval: constant value of the integrand
%   order: spline order
%   ngauss: number of gauss points
%   di,dj: derivative order of the splines
%   const: constraints of the splines
%   xknots_core: spline knots for H-mode spline basis on [0,rhoped]

% F. Felici CRPP/EPFL 2010

% check optional input
nsp = numel(xknots)+order-2;
T = eye(nsp);

xknots_core = [];
if ~isempty(varargin)
  xknots_core=varargin{1};
end

nx = length(xknots);
[xgauss,wgauss] = GaussGrid(xknots,ngauss);

if numel(fval) == 1;
  fval = fval*ones(size(xgauss)); % assume constant fval;
elseif numel(fval)~=numel(xgauss);
  error('fval must be given on the xgauss points')
end

if isempty( xknots_core )
 % Standard case:
	[Lam,Lamp] = splineval(xgauss,xknots,order,constr); % spline values and derivatives at gauss points
else
  % Projection on H mode spline basis on [0,rhoped].
  % Value of linear splines on [rhoped,1] can be set to zero as
  % multiplication with the projection splines would cancel these values
  % anyway.
  [Lam,Lamp] = splineval_hmode(xgauss,xknots_core,order);
  Lam(end-1, xgauss>xknots_core(end)) = 0;
  Lam = Lam(1:end-1, :);
  Lamp(end-1, xgauss>xknots_core(end)) = 0;
  Lamp = Lamp(1:end-1, :);
end

switch di
  case 0
    Lamdi = Lam;
    Lamdj = Lam;
  case 1
    Lamdi = Lamp;
    Lamdj = Lamp;
  otherwise
    error('not implemented yet for di or dj>1')
end

% special case, give it a try
if ~isempty(dj)
  if (di==1 && dj==0)
    Lamdi = Lamp;
    Lamdj = Lam;
  elseif (di==0 && dj==1)
    Lamdi = Lam;
    Lamdj = Lamp;
  end
end

XI = [];
XIindex = [];
if ~isempty(dj) % full FE matrix
  %if di~=dj; error('not implemented yet for di~=dj'); end
  
  if di==dj % symmetric, do upper diagonals only
    for id = 0:order;% diagonal number
      xi = zeros(nsp-id,ngauss*(nx-1));
      xiij = zeros(nsp-id,2);
      for ir = 1:nsp-id
        xi(ir,:) = full([wgauss'.*fval'.*Lamdi(ir,:).*Lamdj(ir+id,:)]);
        % indices of the FE matrix on which this element should end up
        xiij(ir,:) = [ir+id,ir]; % i,j indices
      end
      XI = [XI;full(xi)];
      XIindex = [XIindex;xiij];
    end
    idiag = find(XIindex(:,1) ~= XIindex(:,2)); % indices of off-diagonal terms
    XIindex = [XIindex;fliplr(XIindex(idiag,:))]; % to get both upper and lower diagonal
    XI = [XI;XI(idiag,:)];
  else % asymmetric
    for id = -order:order ;% diagonal number, must do them all
      xi = zeros(nsp-abs(id),ngauss*(nx-1));
      xiij = zeros(nsp-abs(id),2);
      for ir = 1:nsp-abs(id)
        if id>0
          idj = id; idi = 0;
        elseif id<0
          idj = 0; idi = -id;
        else
          idi = 0; idj = 0;
        end
        xi(ir,:) = full([wgauss'.*fval'.*Lamdi(ir+idi,:).*Lamdj(ir+idj,:)]);
        % indices of the FE matrix on which this element should end up
        xiij(ir,:) = [ir+idi,ir+idj]; % i,j indices
      end
      XI = [XI;full(xi)];
      XIindex = [XIindex;xiij];
    end
  end
  
  % linear index:
  XIindexo = sub2ind([nsp,nsp],XIindex(:,1),XIindex(:,2));
  
  % Apply transformation such that
  % (A(XIindex))*T = (XI*p)*T = B(XIindexnew) = XInew*p;
  
  XIn = zeros(nsp.^2,size(XI,2)); % init
  for ix=1:numel(XIindexo)
    %%
    XX=zeros(nsp,nsp);
    % location in original matrix
    XX(XIindexo(ix)) = 1;
    % transformation
    XT = XX*T;
    % nonzero indices after transformation
    ixt = find(XT);
    
    % add old row of XI to this index
    XIn(ixt,:) = bsxfun(@plus,XIn(ixt,:),XI(ix,:));
  end
  
  % store only nonzero entries
  XIindexout = find(any(XIn'));
  XIout = XIn(XIindexout,:);
  
  %% test
  check=0;
  if check
    A = zeros(nsp);
    A(XIindexout) = XIout*ones(size(XIout,2),1);
    subplot(221); imagesc(A)
    
    B = zeros(nsp);
    B(XIindexo) = XI*ones(size(XIout,2),1);
    subplot(222); imagesc(B)
    
    subplot(223); imagesc(B*T)
    subplot(224); imagesc(A-B*T); colorbar
  end
  
  
else % source term only -> vector = XI{1}*source(t)
  xi = zeros(nsp,ngauss*(nx-1));
  % apply spline transformation directly here
  TLamdi = T'*Lamdi;
  for ir = 1:nsp
    xi(ir,:) = full([wgauss'.*fval'.*(TLamdi(ir,:))]);
  end
  XIout = full(xi); % contains all the matrices to construct the diagonals
  XIindexout = [1:nsp]'; % trivial case
end



return
%%



