function [sln,CofQhg] = DarcyGen_WG_QuadriP0P0AC0_AsmSlv( ... 
  EqnBC, BndryCondType, QuadriMesh, PermK, GAUSSQUAD)
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, ColoState; 2014/07--2020/08 

%% Assembly is now separated 
[GlbMat,GlbRHS, GlbVecEss, FlagDOFs] ... 
  = DarcyGen_WG_QuadriP0P0AC0_Asm( ... 
    EqnBC, BndryCondType, QuadriMesh, GAUSSQUAD);

%% Enforcing essential boundary conditions 
sln = GlbVecEss;
GlbRHS = GlbRHS - GlbMat*sln;

%% Solving the reduced global linear system directly 
FreeDOFs = find(~FlagDOFs);
sln(FreeDOFs) = GlbMat(FreeDOFs,FreeDOFs) \ GlbRHS(FreeDOFs);

%% Shown below: Auxiliary but needed stuff 
NumEms = QuadriMesh.NumEms;
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;
GM = Hdiv_QuadriAC0_NmlzPiolaBas_GramMat(QuadriMesh,GAUSSQUAD);

%% Darcy: Computing local L2-projection of g=-(K/mu)f into elementwise AC0 
CofQhg = zeros(NumEms,4);
EltSrsg = zeros(NumEms,4);
MatJ = zeros(NumEms,2,2);
MatP = zeros(NumEms,2,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  J11 = CofA(:,2) + CofA(:,4)*yhat;  
  J12 = CofA(:,3) + CofA(:,4)*xhat;
  J21 = CofB(:,2) + CofB(:,4)*yhat;  
  J22 = CofB(:,3) + CofB(:,4)*xhat;
  jac = J11.*J22 - J12.*J21;
  MatJ(:,1,1) = J11;  MatJ(:,1,2) = J12;
  MatJ(:,2,1) = J21;  MatJ(:,2,2) = J22;
  for i=1:2
    for j=1:2
      MatP(:,i,j) = MatJ(:,i,j)./jac;
    end
  end
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  X = qp(:,1) - QuadriMesh.EmCntr(:,1);
  Y = qp(:,2) - QuadriMesh.EmCntr(:,2);
  U = MatP(:,1,1) .* xhat - MatP(:,1,2) .* yhat;
  V = MatP(:,2,1) .* xhat - MatP(:,2,2) .* yhat;
  gval = EqnBC.fxng(qp);
  EltSrsg(:,1) = EltSrsg(:,1) +  gval(:,1) .* jac *  GAUSSQUAD.RECT(k,3);
  EltSrsg(:,2) = EltSrsg(:,2) +  gval(:,2) .* jac *  GAUSSQUAD.RECT(k,3);
  EltSrsg(:,3) = EltSrsg(:,3) + (gval(:,1).*X + gval(:,2).*Y) .* jac *  GAUSSQUAD.RECT(k,3);
  EltSrsg(:,4) = EltSrsg(:,4) + (gval(:,1).*U + gval(:,2).*V) .* jac *  GAUSSQUAD.RECT(k,3);
end
for ie=1:NumEms 
  CofQhg(ie,:) = squeeze(GM(ie,:,:)) \ squeeze(EltSrsg(ie,:))';
end

return;