function GMK = Hdiv_PlygnCW0_Wachspress_GramMatK(PlygnMesh,MatK,GAUSSQUAD)
%% Hdiv: CW0 elements on plygn.mesh: Weighted Gram matrix GMK 
% Using Wachspress coordinates and normalized coordinates 
% See Chen,Wang, Math. Comput., 2017, p.2060 formula (3.2) 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/08 

%% Mesh info 
NumEms = PlygnMesh.NumEms; 
xc = PlygnMesh.EmCntr(:,1);
yc = PlygnMesh.EmCntr(:,2);

%% Conversion matrix 
% ConvMat = cell(NumEms,1);
ConvMat = Hdiv_PlygnCW0_Wachspress_ConvMatBasFrm(PlygnMesh);

%% Computing weighted Gram matrix by a chosen Gaussian quadrature 
GMK = cell(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  GMK{ie} = zeros(n,n);
  bas = zeros(n,2);
  frm = zeros(n+1,2);
  vrtx = PlygnMesh.node(PlygnMesh.elem{ie}(:),:);
  for m=1:(n-2)
    x1 = PlygnMesh.node(PlygnMesh.elem{ie}(1),1);
    y1 = PlygnMesh.node(PlygnMesh.elem{ie}(1),2);
    x2 = PlygnMesh.node(PlygnMesh.elem{ie}(m+1),1);
    y2 = PlygnMesh.node(PlygnMesh.elem{ie}(m+1),2);
    x3 = PlygnMesh.node(PlygnMesh.elem{ie}(m+2),1);
    y3 = PlygnMesh.node(PlygnMesh.elem{ie}(m+2),2);
    ar = abs(0.5*((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1)));
    for k=1:NumQuadPts
      qp = GAUSSQUAD.TRIG(k,1) * [x1,y1] ...
         + GAUSSQUAD.TRIG(k,2) * [x2,y2] ...
         + GAUSSQUAD.TRIG(k,3) * [x3,y3] ;
      % Normalized coordinates 
      X = qp(1) - xc(ie);
      Y = qp(2) - yc(ie);
      % Wachspress coordinates 
      [lambda,lambdagrad,lambdacurl] = Plygn_Wachspress_ValGradCurl(vrtx,qp);
      % Frame function values 
      frm(1,1) = X;
      frm(1,2) = Y;
      for i=1:n
        frm(1+i,1) = lambdacurl(i,1);
        frm(1+i,2) = lambdacurl(i,2);
      end
      % Basis function values 
      for i=1:n
        bas(i,1) = ConvMat{ie}(i,:)*frm(:,1);
        bas(i,2) = ConvMat{ie}(i,:)*frm(:,2);
      end
      % Integration for GMK 
      for i=1:n
        for j=1:n
          GMK{ie}(i,j) = GMK{ie}(i,j) + ar * GAUSSQUAD.TRIG(k,4) ...
            * ( MatK(ie,1,1) * bas(i,1) * bas(j,1) ...
              + MatK(ie,1,2) * bas(i,1) * bas(j,2) ...
              + MatK(ie,2,1) * bas(i,2) * bas(j,1) ...
              + MatK(ie,2,2) * bas(i,2) * bas(j,2));
        end
      end
    end
  end
end

return;