function [CofDP1,CofRT0] = Darcy_DG_TriP1_ProjPresVel(...
  EqnBC,TriMesh,GAUSSQUAD)
%% Darcy: DG: TriP1: Projecting exact pressure & velocity (assuming known) 
%   CofDP1: pressure -> DP1   CofRT0: velocity -> RT0(eg.bas.,local) 
% James Liu, ColoState; 2012/07--2016/05 

%% Mesh info 
% Basic mesh info 
NumEms = TriMesh.NumEms;
% More mesh info 
area = TriMesh.area;
% k1 = TriMesh.elem(:,1);  k2 = TriMesh.elem(:,2);  k3 = TriMesh.elem(:,3);
% x1 = TriMesh.node(k1,1);  y1 = TriMesh.node(k1,2); 
% x2 = TriMesh.node(k2,1);  y2 = TriMesh.node(k2,2);
% x3 = TriMesh.node(k3,1);  y3 = TriMesh.node(k3,2);

% %% Computing emws. Gram matrix for DP1: Approach I 
% IX2 = (area/36).*((x1-x2).^2+(x2-x3).^2+(x3-x1).^2);
% IY2 = (area/36).*((y1-y2).^2+(y2-y3).^2+(y3-y1).^2);
% % JL20131101: Formula for IXY might be incorrect 
% IXY = (area/36).*(2*(x2-x1).*(y2-y1)-(x2-x1).*(y3-y1)...
%                  -(x3-x1).*(y2-y1)+2*(x3-x1).*(y3-y1));
% GMDP1 = zeros(NumEms,3,3);  % emws. Gram matrices for DP1 
% GMDP1(:,1,1) = area;  GMDP1(:,2,2) = IX2;  GMDP1(:,3,3) = IY2;
% GMDP1(:,2,3) = IXY;  GMDP1(:,3,2) = IXY;

%% Computing emws. Gram matrix for DP1: Approach II 
GMDP1 = zeros(NumEms,3,3);  % emws. Gram matrices for DP1 
WO = ones(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1)*TriMesh.node(TriMesh.elem(:,1),:)...
     + GAUSSQUAD.TRIG(k,2)*TriMesh.node(TriMesh.elem(:,2),:)...
     + GAUSSQUAD.TRIG(k,3)*TriMesh.node(TriMesh.elem(:,3),:);
  X = qp(:,1);  Y = qp(:,2);
  BasFxn = [WO,X,Y];
  for i=1:3
    for j=i:3
      GMDP1(:,i,j) = GMDP1(:,i,j)...
        + GAUSSQUAD.TRIG(k,4)*(BasFxn(:,i).*BasFxn(:,j)).*area;
      GMDP1(:,j,i) = GMDP1(:,i,j);  % Utilizing symmetry 
    end
  end
end

%% Projecting exact pressure to discontinuous poly. deg.1 (DP1) 
% Computing emws. right-hand side for DP1 
RHS = zeros(NumEms,3);
WO = ones(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1)*TriMesh.node(TriMesh.elem(:,1),:)...
     + GAUSSQUAD.TRIG(k,2)*TriMesh.node(TriMesh.elem(:,2),:)...
     + GAUSSQUAD.TRIG(k,3)*TriMesh.node(TriMesh.elem(:,3),:);
  X = qp(:,1);  Y = qp(:,2);
  BasFxn = [WO,X,Y];
  pval = EqnBC.fxnp(qp);
  for j=1:3
    RHS(:,j) = RHS(:,j) + GAUSSQUAD.TRIG(k,4)*(pval.*BasFxn(:,j)).*area;
  end
end
% Solving lin.sys. for coeff. of projected pressure in DP1 natural basis 
CofDP1 = zeros(NumEms,3);
for ie=1:NumEms
  EltGM = squeeze(GMDP1(ie,:,:));
  EltRHS = RHS(ie,:)';
  CofDP1(ie,:) = (EltGM\EltRHS)';
end

%% Projecting exact velocity to RT0 (edge based basis functions, local) 
% Gram matrices of RT0 eg.bas.bas.fxns.(local) 
GMRT0 = Hdiv_TriRT0_GramMat_EgBas(TriMesh);
% Computing "RHS" 
RHS = zeros(NumEms,3);
BasFxn = zeros(NumEms,2,3);
NumQuadPts = size(GAUSSQUAD.TRIG,1); 
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1)*TriMesh.node(TriMesh.elem(:,1),:)...
     + GAUSSQUAD.TRIG(k,2)*TriMesh.node(TriMesh.elem(:,2),:)...
     + GAUSSQUAD.TRIG(k,3)*TriMesh.node(TriMesh.elem(:,3),:);
  vel = -repmat(EqnBC.fxnK(qp),1,2).*EqnBC.fxnpg(qp);
  for j=1:3
    coeff = TriMesh.LenEg(TriMesh.elem2edge(:,j))./(2*area);
    BasFxn(:,:,j) = [coeff,coeff].*(qp-TriMesh.node(TriMesh.elem(:,j),:));
    RHS(:,j) = RHS(:,j)...
      + GAUSSQUAD.TRIG(k,4)*dot(vel,BasFxn(:,:,j),2).*area;
  end
end
% Solving lin.sys. for coeff. of proj.vel. in RT0 eg.bas.bas.fxns.(local)
CofRT0 = zeros(NumEms,3);  % all elements, 3 edges, 2 ends 
for ie=1:NumEms
  EltGM = squeeze(GMRT0(ie,1:3,1:3));
  EltRHS = RHS(ie,1:3)';
  CofRT0(ie,:) = (EltGM\EltRHS)';
end

return;