Globals
r = JacobiGL(0,0,N);
V  = Vandermonde1D(N, r); invV = inv(V);

tempLIFT = Lift1D; LIFT=[];  for k=1:K; LIFT=blkdiag(LIFT,tempLIFT); end;
Vr = zeros(length(r),(N+1));for i=0:N;   [Vr(:,i+1)] = GradJacobiP(r(:),0,0,i); end;
tempDr = Vr/V; Dr=[]; for k=1:K; Dr=blkdiag(Dr,tempDr); end;
va = EToV(:,1)'; vb = EToV(:,2)'; z = ones(N+1,1)*VZ(va) + 0.5*(r+1)*(VZ(vb)-VZ(va));

zr  = tempDr*z; rz = 1./zr;
Dz = Dr*diag(rz(:)); % global DG differentiation operator
fmask1 = find( abs(r+1) < 1e-10)'; fmask2 = find( abs(r-1) < 1e-10)'; 
Fmask  = [fmask1;fmask2]';
nz=[]; nz(1,1)=-1; nz(2,1)=1; nz = nz*ones(1,K); nz=nz(:); 
Fscale = 1./(zr(Fmask,:)); Fscale = Fscale(:);
vmapM = zeros(2*K,1); vmapP = zeros(2*K,1);
nodeids = reshape(1:Np*K, K*Np,1);
for k=1:K; 
    aaa=nodeids((k-1)*Np+1:k*Np,:); vmapM(2*(k-1)+1:2*k,:) = (aaa(Fmask,:)); 
end;
for k=1:K; 
    km=vmapM(2*(k-1)+1,:); kp=vmapM(2*k,:); 
    vmapP(2*(k-1)+1,:) = km-1; vmapP(2*k,:) = kp+1; 
end;
vmapP(1,:) = vmapM(1,:); vmapP(end,:) = vmapM(end,:);