%% JL20131101: TO BE FINISHED 
function [CofDP2,CofBDM1] = Darcy_DG_TriP2_ProjPresVel(...
  EqnBC,TriMesh,GAUSSQUAD)
%% Darcy: DG: TriP2: Projecting exact pressure & velocity (assuming known) 
%   CofP2: pressure -> DP2   CofBDM1: velocity -> BDM1(local)
% James Liu, ColoState; 2012/07--2016/05 

%% Mesh info 
% Basic mesh info 
NumEms = TriMesh.NumEms;
% Assuming TriMesh.flag>=2 
area = TriMesh.area;

%% Projecting exact pressure to DP2 
% Computing emws. Gram matrices and right-hand sides for DP2 
GMDP2 = zeros(NumEms,6,6);  % emws. Gram matrices for DP2 
RHS = zeros(NumEms,6);
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, X.^2, X.*Y, Y.^2];
  for i=1:6
    for j=i:6
      GMDP2(:,i,j) = GMDP2(:,i,j)...
        + GAUSSQUAD.TRIG(k,4)*(BasFxn(:,i).*BasFxn(:,j)).*area;
      GMDP2(:,j,i) = GMDP2(:,i,j);  % Utilizing symmetry 
    end
  end
  pval = EqnBC.fxnp(qp);
  for j=1:6
    RHS(:,j) = RHS(:,j) + GAUSSQUAD.TRIG(k,4)*(pval.*BasFxn(:,j)).*area;
  end
end
% Solving lin.sys. for coeff. of projected pressure in natural basis 
CofDP2 = zeros(NumEms,6);
for ie=1:NumEms
  EltGM = squeeze(GMDP2(ie,:,:));
  EltRHS = RHS(ie,:)';
  CofDP2(ie,:) = (EltGM\EltRHS)';
end

%% Auxiliary quantities 
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);
eg1 = [x3,y3] - [x2,y2];
eg2 = [x1,y1] - [x3,y3];
eg3 = [x2,y2] - [x1,y1];
lg1 = sqrt(eg1(:,1).^2+eg1(:,2).^2);
lg2 = sqrt(eg2(:,1).^2+eg2(:,2).^2);
lg3 = sqrt(eg3(:,1).^2+eg3(:,2).^2);
tan = zeros(NumEms,2,3);
tan(:,:,1) = eg1./[lg1,lg1];
tan(:,:,2) = eg2./[lg2,lg2];
tan(:,:,3) = eg3./[lg3,lg3];
% cosine theta 
ct1 = (lg2.^2+lg3.^2-lg1.^2)./(2*lg2.*lg3);
ct2 = (lg3.^2+lg1.^2-lg2.^2)./(2*lg3.*lg1);
ct3 = (lg1.^2+lg2.^2-lg3.^2)./(2*lg1.*lg2);
% sine theta 
st1 = sqrt(1-ct1.^2);
st2 = sqrt(1-ct2.^2);
st3 = sqrt(1-ct3.^2);

%% Projecting exact velocity to BDM1 
% Gram matrices of BDM1 eg.bas.bas.fxns.(local) 
GMBDM1 = Hdiv_TriBDM1_GramMat_EgBas(TriMesh);
% Computing "RHS" 
RHS = zeros(NumEms,6);
BasFxn = zeros(NumEms,2,6);
NumQuadPts = size(GAUSSQUAD.TRIG,1); 
% JL20131017: TO BE FINISHED 
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),:);
  pg = EqnBC.fxnpg(qp);
  K = EqnBC.fxnK(qp);
  vel = -[K,K].*pg;
  BaryCrd = TriMesh_BaryCrd(TriMesh,qp);
  BasFxn(:,:,1) =  repmat(BaryCrd(:,2)./st2,1,2).*tan(:,:,3);
  BasFxn(:,:,2) = -repmat(BaryCrd(:,3)./st3,1,2).*tan(:,:,2);
  BasFxn(:,:,3) =  repmat(BaryCrd(:,3)./st3,1,2).*tan(:,:,1);
  BasFxn(:,:,4) = -repmat(BaryCrd(:,1)./st1,1,2).*tan(:,:,3);
  BasFxn(:,:,5) =  repmat(BaryCrd(:,1)./st1,1,2).*tan(:,:,2);
  BasFxn(:,:,6) = -repmat(BaryCrd(:,2)./st2,1,2).*tan(:,:,1);
  for j=1:6
    ww = squeeze(BasFxn(:,:,j));
    RHS(:,j) = RHS(:,j) + GAUSSQUAD.TRIG(k,4)*dot(vel,ww,2).*area;
  end
end
% Solving lin.sys. for coeffs. of proj.vel. in BDM1 eg.bas.bas.fxns.
CofBDM1 = zeros(NumEms,6);  % all elements, 3 edges, 2 ends 
for ie=1:NumEms
  EltGM = squeeze(GMBDM1(ie,:,:));
  EltRHS = RHS(ie,:)';
  CofBDM1(ie,1:6) = (EltGM\EltRHS)';
end

return;