function [val] = FEval(xvec,xgrid,i,derorder)
%function [val] = FEval(xvec,xgrid,i,derorder)
% value for pwl functions
% xvec vector to evaluate value
% xgrid grid of FE nodes
% i element required
% derorder = 0 for value, derorder = 1 to return derivative

    if any(diff(xgrid)<=0); error('xgrid not increasing'); end

    xn = xgrid;
    mm1 = ones(size(xvec));
    if i==1 % first node
        hM = mm1.*(xn(i+1)-xn(i)); 
        lM = [(mm1.*xn(i+1)-xvec)];
        signX = 0.5*sign(xvec-xn(i))+0.5;
        hi = (signX.*hM);
        li = signX.*lM;
    elseif i==length(xn); % last node    
        hM =  mm1.*(xn(i)-xn(i-1)); 
        lM = [(xvec-mm1.*xn(i-1))];
        signX = 0.5*sign(xn(i)-xvec)+0.5;
        hi = (signX.*hM);
        li = signX.*lM;
    else
        hM = [(mm1.*xn(i+1)-xn(i));...
            (mm1.*xn(i)-xn(i-1))]; 
        lM = [(mm1.*xn(i+1)-xvec);...
            (xvec-mm1.*xn(i-1))];
        signX = [0.5*sign(xvec-xn(i))+0.5;
            0.5*sign(xn(i)-xvec)+0.5];

        hi = ([1 1]*(signX.*hM));
        li = ([1 1]*(signX.*lM));

    end
    xi = (xvec-xn(i))./hi;
    L1 = li./hi; % piecewise linear
    L1p = -sign(xi)./hi; % derivative
    
    G0  = (abs(xi) - 1).^2.*(2*abs(xi)+1); % first hermite
    G1  = (abs(xi) - 1).^2.*xi.*hi; % second hermite
    G0p = (6*abs(xi).^2.*sign(xi)-6*abs(xi).*sign(xi))./hi;
    G1p = (abs(xi)-1).*(abs(xi)+2*xi.*sign(xi)-1);
    
    jthisel = find(abs(xi)<=1);
    
    L1(abs(xi)>1) = 0;
    L1p(abs(xi)>1) = 0;
    
    if derorder == 0
        val = L1;
    elseif derorder == 1;
        val = L1p;
    else
        error('invalid derorder')
    end

return