function [sln,RedGlbMat] = Poisson_WG_QuadriP0P0AC0Schur_AsmSlv(...
  EqnBC,BndryCondType,QuadriMesh,CofK,GAUSSQUAD)
%% Poisson: WG(P0,P0;AC0) on a quadrilateral mesh: Assembly & Solving 
% Using local Schur complement 
% Assuming CofK is an elementwise constant 2x2 SPD matrix 
% James Liu, ColoState; 2012/07--2021/04 

%% Mesh info 
NumEms = QuadriMesh.NumEms;
NumEgs = QuadriMesh.NumEgs;
% Assuming QuadriMesh.flag>=2 
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% Sorting out boundary edges: Dirichlet, Neumann 
DirichletEdge = find(BndryCondType(QuadriMesh.BndryEdge+1)==1);
NeumannEdge   = find(BndryCondType(QuadriMesh.BndryEdge+1)==2);

%% Boundary info
NumDirichletEgs = size(DirichletEdge,1);
NumNeumannEgs = size(NeumannEdge,1);

%% Auxiliary arrays 
% GMK = zeros(NumEms,4,4);
GMK = Hdiv_QuadriAC0_NmlzPiolaBas_GramMatK(QuadriMesh,CofK,GAUSSQUAD);
% CDWGB = zeros(NumEms,5,4);
CDWGB = WG_QuadriP0P0AC0_CofAC0NmlzPiolaBas_DiscWkGradBasFxn(QuadriMesh,GAUSSQUAD);
EltMat = zeros(NumEms,5,5);
%% JL20171011: TO BE REVISED FOR EFFECIENCY !!
for ie=1:NumEms
  EGM = squeeze(GMK(ie,:,:));
  for i=1:5
    for j=1:5
      EltMat(ie,i,j) = squeeze(CDWGB(ie,i,:))' * EGM * squeeze(CDWGB(ie,j,:));
    end
  end
end

%% Performing local/elementwise Schur complement 
A00 = squeeze(EltMat(:,1,1));
A01 = squeeze(EltMat(:, 1, 2:5));
A10 = squeeze(EltMat(:, 2:5, 1));
A11 = EltMat(:, 2:5, 2:5);
S11 = zeros(NumEms,4,4);
for i=1:4
  for j=1:4
    S11(:,i,j) = A11(:,i,j) - A10(:,i) ./ A00(:,1) .* A01(:,j);
  end
end

%% Setting up the global linear system for the mesh skeleton 
DOFs = NumEgs;
% GlbMat = sparse(DOFs,DOFs);
% GlbRHS = zeros(DOFs,1);

%% Assembling the global coefficient matrix GlbMat 
% Positions 
pos = zeros(NumEms,4);
for j=1:4
  pos(:,j) = QuadriMesh.elem2edge(:,j);
end
% Now assembly 
GlbMat = sparse(DOFs,DOFs);
for i=1:4
  II = pos(:,i); 
  for j=1:4
    JJ = pos(:,j);
    GlbMat = GlbMat + sparse(II,JJ,S11(:,i,j),DOFs,DOFs);
  end
end

%% Assembling GlbRHS 
GlbRHS = zeros(DOFs,1);
% NumQuadPts = size(GAUSSQUAD.RECT,1);
% GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
% qp = zeros(NumEms,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);
%   GlbRHS(1:NumEms) = GlbRHS(1:NumEms) ...
%     + GAUSSQUAD.RECT(k,3) * jac .* EqnBC.fxnf(qp);
% end

% %% Incorporating boundary conditions: Neumann as natural 
% if NumNeumannEgs>0
%   x1 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,1),1);
%   y1 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,1),2);
%   x2 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,2),1);
%   y2 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,2),2);
%   LenNeumannEg = sqrt((x2-x1).^2+(y2-y1).^2);
%   NumQuadPts = size(GAUSSQUAD.LINE,1);
%   qp = zeros(NumNeumannEgs,2);
%   for k=1:NumQuadPts
%     qp(:,1) = GAUSSQUAD.LINE(k,1)*x1 + GAUSSQUAD.LINE(k,2)*x2;
%     qp(:,2) = GAUSSQUAD.LINE(k,1)*y1 + GAUSSQUAD.LINE(k,2)*y2;
%     GlbRHS(NumEms+NeumannEdge) = GlbRHS(NumEms+NeumannEdge)...
%       - GAUSSQUAD.LINE(k,3) * EqnBC.fxnuN(qp) .* LenNeumannEg;
%   end
% end

%% For reducing... 
flag = zeros(DOFs,1);
flag(DirichletEdge) = ones(NumDirichletEgs,1);
FreeEg = find(~flag);
sln = zeros(DOFs,1);

%% Incorporating boundary conditions: Dirichlet as essential 
if NumDirichletEgs>0 
  x1 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,1),1);
  y1 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,1),2);
  x2 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,2),1);
  y2 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,2),2);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  qp = zeros(NumDirichletEgs,2);
  pDavg = zeros(NumDirichletEgs,1);
  for k=1:NumQuadPts
    qp(:,1) = GAUSSQUAD.LINE(k,1)*x1 + GAUSSQUAD.LINE(k,2)*x2;
    qp(:,2) = GAUSSQUAD.LINE(k,1)*y1 + GAUSSQUAD.LINE(k,2)*y2;
    pDavg = pDavg + GAUSSQUAD.LINE(k,3) * EqnBC.fxnpD(qp);
  end
  sln(DirichletEdge) = pDavg;
end

%% Reducing...
GlbRHS = GlbRHS - GlbMat*sln;

%% Solving the reduced global linear system directly 
sln(FreeEg) = GlbMat(FreeEg,FreeEg) \ GlbRHS(FreeEg);
RedGlbMat = GlbMat(FreeEg,FreeEg);

return;