function [ProjCofQ,ProjCofT] = Hdiv_HybridRT0_NmlzBas_ProjMatK(...
  HybridMesh,PermKt,GAUSSQUAD)
% Projection coefficients of (MatK * nmlz.bas.fxn.) back in nmlz.bas.fxn. 
% For a hybrid mesh: Quadrilaterals, Triangles 
% MatK is a constant 2x2 matrix on each element 
% ProjCof has size #Elements*4*4 or #Elements*3*3 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/02 

%% Mesh info 
% NumEms = HybridMesh.NumEms;
NumEmsQ = HybridMesh.NumEmsQ;
NumEmsT = HybridMesh.NumEmsT;
LblQ = HybridMesh.labelemQ;
LblT = HybridMesh.labelemT;

%% Auxiliary quantities for quadrilateral elements
k1 = HybridMesh.elemQ(:,1);  k2 = HybridMesh.elemQ(:,2);  
k3 = HybridMesh.elemQ(:,3);  k4 = HybridMesh.elemQ(:,4);
x = zeros(NumEmsQ,5);  y = zeros(NumEmsQ,5);
x1 = HybridMesh.node(k1,1);  y1 = HybridMesh.node(k1,2);
x2 = HybridMesh.node(k2,1);  y2 = HybridMesh.node(k2,2);
x3 = HybridMesh.node(k3,1);  y3 = HybridMesh.node(k3,2);
x4 = HybridMesh.node(k4,1);  y4 = HybridMesh.node(k4,2);
x(:,1) = HybridMesh.node(k1,1);  y(:,1) = HybridMesh.node(k1,2);
x(:,2) = HybridMesh.node(k2,1);  y(:,2) = HybridMesh.node(k2,2);
x(:,3) = HybridMesh.node(k3,1);  y(:,3) = HybridMesh.node(k3,2);
x(:,4) = HybridMesh.node(k4,1);  y(:,4) = HybridMesh.node(k4,2);
x(:,5) = x(:,1);                 y(:,5) = y(:,1);
xc = 0.25 * sum(x(:,1:4),2);
yc = 0.25 * sum(y(:,1:4),2);
% Coefficients for the bilinear mapping for quadrilaterals
CofA = zeros(NumEmsQ,4);      CofB = zeros(NumEmsQ,4);
CofA(:,1) = x1;               CofB(:,1) = y1;
CofA(:,2) = x2-x1;            CofB(:,2) = y2-y1;
CofA(:,3) = x4-x1;            CofB(:,3) = y4-y1;
CofA(:,4) = (x1+x3)-(x2+x4);  CofB(:,4) = (y1+y3)-(y2+y4);

%% Computing integrals for quadrilaterals 
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
S1 = zeros(NumEmsQ,1);  SX2 = zeros(NumEmsQ,1);
SX = zeros(NumEmsQ,1);  SY2 = zeros(NumEmsQ,1);
SY = zeros(NumEmsQ,1);  SXY = zeros(NumEmsQ,1);
qp = zeros(NumEmsQ,2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp(:,1) = CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat;
  qp(:,2) = CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat;
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  X = qp(:,1) - xc;
  Y = qp(:,2) - yc;
  S1 = S1 + GAUSSQUAD.RECT(k,3)*jac;  % S1=area
  SX = SX + GAUSSQUAD.RECT(k,3)*jac.*X;
  SY = SY + GAUSSQUAD.RECT(k,3)*jac.*Y;
  SX2 = SX2 + GAUSSQUAD.RECT(k,3)*jac.*(X.*X);
  SY2 = SY2 + GAUSSQUAD.RECT(k,3)*jac.*(Y.*Y);
  SXY = SXY + GAUSSQUAD.RECT(k,3)*jac.*(X.*Y);
end

%% Gram matrix for quadrilaterlas RT[0] 
GM = zeros(NumEmsQ,4,4);
GM(:,1,1) = S1;  GM(:,2,2) = S1;  % S1=area
GM(:,1,3) = SX;  GM(:,2,4) = SY;
GM(:,3,1) = SX;  GM(:,4,2) = SY;
GM(:,3,3) = SX2;  GM(:,4,4) = SY2;
% Cofs. in quadri.RT[0]nmlz.bas. for projection of permeability * bas.fxn.
ProjCofQ = zeros(NumEmsQ,4,4);  % i-th fxn. projected onto k-th fxn. 
% JL20170116: TO BE REVISED FOR EFFECIENCY! 
for ie=1:NumEmsQ
  EltGM = squeeze(GM(ie,:,:));
  % For the 1st bas.fxn.
  EltRHS = [PermKt(ie,1,1)*S1(ie);...
            PermKt(ie,2,1)*S1(ie);...
            PermKt(ie,1,1)*SX(ie);...
            PermKt(ie,2,1)*SY(ie)];
  ProjCofQ(ie,1,:) = EltGM\EltRHS;
  % For the 2nd bas.fxn.
  EltRHS = [PermKt(ie,1,2)*S1(ie);...
            PermKt(ie,2,2)*S1(ie);...
            PermKt(ie,1,2)*SX(ie);...
            PermKt(ie,2,2)*SY(ie)];
  ProjCofQ(ie,2,:) = EltGM\EltRHS;
  % For the 3rd bas.fxn.
  EltRHS = [PermKt(ie,1,1)*SX(ie);...
            PermKt(ie,2,1)*SX(ie);...
            PermKt(ie,1,1)*SX2(ie);...
            PermKt(ie,2,1)*SXY(ie)];
  ProjCofQ(ie,3,:) = EltGM\EltRHS;
  % For the 4th bas.fxn.
  EltRHS = [PermKt(ie,1,2)*SY(ie);...
            PermKt(ie,2,2)*SY(ie);...
            PermKt(ie,1,2)*SXY(ie);...
            PermKt(ie,2,2)*SY2(ie)];
  ProjCofQ(ie,4,:) = EltGM\EltRHS;
end

%% Mesh info for Triangular elems.
k1 = HybridMesh.elemT(:,1);  k2 = HybridMesh.elemT(:,2);  k3 = HybridMesh.elemT(:,3);
x1 = HybridMesh.elemT(k1,1);  y1 = HybridMesh.elemT(k1,2); 
x2 = HybridMesh.elemT(k2,1);  y2 = HybridMesh.elemT(k2,2);
x3 = HybridMesh.elemT(k3,1);  y3 = HybridMesh.elemT(k3,2);
areaT = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));

%% Auxiliary quantities 
SX2 = (areaT/36).*( (x1-x2).^2 + (x2-x3).^2 + (x3-x1).^2 );
SY2 = (areaT/36).*( (y1-y2).^2 + (y2-y3).^2 + (y3-y1).^2 );
SXY = (areaT/36).*((x1-x2).*(y1-y2) + (x2-x3).*(y2-y3) + (x3-x1).*(y3-y1));
C3 = ( PermKt(LblT,1,1).*SX2 + PermKt(LblT,1,2).*SXY ... 
     + PermKt(LblT,2,1).*SXY + PermKt(LblT,2,2).*SY2)./(SX2+SY2);

%% Projection coefficients 
ProjCofT = zeros(NumEmsT,3,3);
ZO = zeros(NumEmsT,1);
ProjCofT(:,1,:) = [PermKt(LblT,1,1), PermKt(LblT,2,1), ZO];
ProjCofT(:,2,:) = [PermKt(LblT,1,2), PermKt(LblT,2,2), ZO];
ProjCofT(:,3,:) = [ZO, ZO, C3];

return