function [GM] = Hdiv_PlygnCW0_Wachspress_GramMat(PlygnMesh,GAUSSQUAD)
%% Hdiv: CW0 elements on plygn.mesh: Gram matrix: Using Wachspress crd.
% 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);

%% For frame: "Gram matrix" 
GMF = cell(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1); 
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  vrtx = PlygnMesh.node(PlygnMesh.elem{ie}(:),:);
  GMF{ie} = zeros(n+1,n+1);
  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] ;
      X = qp(1) - xc(ie);
      Y = qp(2) - yc(ie);
      [lambda,lambdagrad,lambdacurl] = Plygn_Wachspress_ValGradCurl(vrtx,qp);
      GMF{ie}(1,1) = GMF{ie}(1,1) + GAUSSQUAD.TRIG(k,4) * (X.*X+Y.*Y) * ar;
      for i=1:n
        GMF{ie}(1+i,1) = GMF{ie}(1+i,1) ...
          + GAUSSQUAD.TRIG(k,4) * (lambdacurl(i,1)*X+lambdacurl(i,2)*Y) * ar;
        GMF{ie}(1,1+i) = GMF{ie}(1+i,1);
        for j=1:n
          GMF{ie}(1+i,1+j) = GMF{ie}(1+i,1+j) ...
            + GAUSSQUAD.TRIG(k,4) * (lambdacurl(i,1)*lambdacurl(j,1) + lambdacurl(i,2)*lambdacurl(j,2)) * ar;
        end
      end
    end
  end
end

%%
GM = cell(NumEms,1);
for ie=1:NumEms
  C = ConvMat{ie};
  GM{ie} = C * GMF{ie} * C';
end

return;