function [L2ErrPres,L2ErrVel,L2ErrFlux,L2ErrDiv] = ...
  Darcy_WG_PlygnP0P0CW0_Err( ...
  EqnBC,PlygnMesh,NumerPresCofP0,NumerVelCofCW0,GAUSSQUAD)
%% Darcy: WG(P0,P0;CW0)Plygn: Computing errors 
% James Liu, ColoState; 2012/07--2018/12 

%% Mesh info 
NumEms = PlygnMesh.NumEms;

%% All kinds of error 
%     L2ErrPres: L^2-norm of error in pressure (for element interior) 
%      L2ErrVel: L^2-norm of error in velocity 
%     L2ErrFlux: L^2-norm of error in flux 
%      L2ErrDiv: L^2-norm of error in divergence 

%% Computing L2-norm of error in pressure (for element interiors) 
ErrPresEm = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  ErrPresEm(ie) = 0;
  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] ;
      pval = EqnBC.fxnp(qp);
      perr = pval - NumerPresCofP0(ie);
      ErrPresEm(ie) = ErrPresEm(ie) + GAUSSQUAD.TRIG(k,4)*(perr^2)*ar;
    end
  end
end
L2ErrPres = sqrt(sum(ErrPresEm));

%% Auxiliary  
ConvMat = Hdiv_PlygnCW0_Wachspress_ConvMatBasFrm(PlygnMesh);

%% Computing L2-norm of error in velocity 
% Using NumerVelCofCW0 
ErrVelEm = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  vrtx = PlygnMesh.node(PlygnMesh.elem{ie}(:),:);
  ErrVelEm(ie) = 0;
  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] ;
      uval = EqnBC.fxnu(qp);
      X = qp(1) - PlygnMesh.EmCntr(ie,1);
      Y = qp(2) - PlygnMesh.EmCntr(ie,2);
      [lambda,lambdagrad,lambdacurl] = ... 
        Plygn_Wachspress_ValGradCurl(vrtx,qp);
      frm = zeros(n+1,2);
      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
      bas = zeros(n,2);
      for i=1:n
        bas(i,1) = ConvMat{ie}(i,:) * frm(:,1);
        bas(i,2) = ConvMat{ie}(i,:) * frm(:,2);
      end
      NV = zeros(1,2);
      NV(1) = NumerVelCofCW0{ie} * bas(:,1);
      NV(2) = NumerVelCofCW0{ie} * bas(:,2);
      uerr = uval - NV;
      ErrVelEm(ie) = ErrVelEm(ie) ... 
        + GAUSSQUAD.TRIG(k,4) * (uerr(1)^2+uerr(2)^2) * ar;
    end
  end
end
L2ErrVel = sqrt(sum(ErrVelEm));

%% Computing L2-norm of error in flux 
ErrFlux = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.LINE,1);
for ie=1:NumEms
  ErrFlux(ie) = 0;
  n = PlygnMesh.ElemType(ie);
  for m=1:n
    m1 = mod(m,n) + 1;
    x1 = PlygnMesh.node(PlygnMesh.elem{ie}(m),1);
    y1 = PlygnMesh.node(PlygnMesh.elem{ie}(m),2);
    x2 = PlygnMesh.node(PlygnMesh.elem{ie}(m1),1);
    y2 = PlygnMesh.node(PlygnMesh.elem{ie}(m1),2);
    PseudoNml = [y2-y1,x1-x2];  % Normal vector carrying with edge length 
    LenEg = sqrt(PseudoNml(:,1)^2+PseudoNml(:,2)^2);
    nml = PseudoNml./[LenEg,LenEg];  % Now unit normal vector 
    NumerFlux = NumerVelCofCW0{ie}(m);
    for k=1:NumQuadPts
      qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
      uval = EqnBC.fxnu(qp);
      ExactFlux = sum(uval.*nml,2);
      ferr = ExactFlux - NumerFlux;
      ErrFlux(ie) = ErrFlux(ie) + (ferr^2)*GAUSSQUAD.LINE(k,3);  % */LenEg;
    end
  end
  ErrFlux(ie) = ErrFlux(ie) * PlygnMesh.area(ie);
end
L2ErrFlux = sqrt(sum(ErrFlux));

%% Computing L2-norm of error in divergence 
% Note: Exact divergence = source: div(u) = f in Darcy equation 
ErrDivEm = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  NumerDiv = NumerVelCofCW0{ie} * ConvMat{ie}(:,1) * 2;  % 2 b/c nmlz.crd.
  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] ;
      divu = EqnBC.fxnf(qp);
      derr = divu - NumerDiv;
      ErrDivEm(ie) = ErrDivEm(ie) + GAUSSQUAD.TRIG(k,4) * (derr.^2) * ar;
    end
  end
end
L2ErrDiv = sqrt(sum(ErrDivEm));

return;
