function Nk = FEMinterpol(X,ngrad)
  if nargin==1,
    ngrad=0;
  end;
  
  nk{1} =inline('(0.5*(x.*x-x))');
  nk{2} =inline('(1-x.*x)');
  nk{3} =inline('(0.5*(x.*x+x))');
  
  dnk{1} = inline('(x-0.5)');
  dnk{2} = inline('((-2)*x)');
  dnk{3} = inline('(x+0.5)');
  
  ix = Q27index+2;   %call the holy function that depends on point order
                     %on a 27 Q quadratic element.
  
  
  % let say I have a column of locations (xi,eta,zeta) in X, then I compute
  % the matrix Nk(k,l), where k runs over the points (along a column) and l
  % runs horizontally along the index of the interpolation function. 
  % in other words, Nk(k,n) contains the n-th function evaluated at the
  % k-th domain point
  
  xi = X(:,1);
  eta= X(:,2);
  zeta=X(:,3);
  
  switch ngrad,
   case 0,
    for k=1:27,
      Nk(:,k) = nk{ix(k,1)}(xi) .*  nk{ix(k,2)}(eta) .* nk{ix(k,3)}(zeta);
    end;
   case 1,
    for k=1:27,
      Nk(:,k) = dnk{ix(k,1)}(xi) .*  nk{ix(k,2)}(eta) .* nk{ix(k,3)}(zeta);
    end;
   case 2,
    for k=1:27,
      Nk(:,k) = nk{ix(k,1)}(xi) .*  dnk{ix(k,2)}(eta) .* nk{ix(k,3)}(zeta);
    end;
   case 3,
    for k=1:27,
      Nk(:,k) = nk{ix(k,1)}(xi) .*  nk{ix(k,2)}(eta) .* dnk{ix(k,3)}(zeta);
    end;
    
   otherwise,
    disp('wrong index in FEMinterpol');
    Nk=[];
  end;