function [sln] = Darcy_WG_PlygnP0P0CW0_AsmSlv( ...
  EqnBC,BndryCondType,PlygnMesh,PermK,GAUSSQUAD)
%% Darcy: WG(P0,P0;CW0) on a polygonal mesh: Assembly & Solving 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2018/12 

%% Mesh info 
NumEms = PlygnMesh.NumEms;
NumEgs = PlygnMesh.NumEgs;

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

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

%% Pre-assembly 
% GMK cell NumEms*1; Each entry n*n SPD matrix; n=ElemType varies 
GMK = Hdiv_PlygnCW0_Wachspress_GramMatK(PlygnMesh,PermK,GAUSSQUAD);
% CDWGB cell NumEms*1; Each entry (n+1)*n matrix; n=ElemType varies 
CDWGB = WG_PlygnP0P0CW0_Wachspress_DiscWkGradBasFxn(PlygnMesh,GAUSSQUAD);

%% NOTE: FACT: Polygon #sides<=maxn 
maxn = PlygnMesh.maxn;
%% Elementwise arrays 
ArrayEE = zeros(NumEms,1);
ArrayEG = zeros(NumEms,maxn);
ArrayGG = zeros(NumEms,maxn,maxn);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  EGGMK = GMK{ie};
  ArrayEE(ie) = CDWGB{ie}(1,:) * EGGMK * (CDWGB{ie}(1,:))';
  for i=1:n
    ArrayEG(ie,i) = CDWGB{ie}(1,:) * EGGMK * (CDWGB{ie}(1+i,:))';
    for m=1:n
      ArrayGG(ie,i,m) = CDWGB{ie}(1+i,:) * EGGMK * (CDWGB{ie}(1+m,:))';
    end
  end
end

%% Setting up the global linear system 
DOFs = NumEms + NumEgs;
% GlbMat = sparse(DOFs,DOFs);
% GlbRHS = zeros(DOFs,1);

%% NOTE: FACT: Polygon #sides<=maxn 
%% Special treatment for edge position (gpos) for assembly 
gpos = DOFs * ones(NumEms,maxn);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  for m=1:n
    gpos(ie,m) = NumEms + PlygnMesh.elem2edge{ie}(m);
  end
end

%% NOTE: FACT: Polygon #sides<=maxn 
%% Assembling the global coefficient matrix
GlbMat = sparse(DOFs,DOFs);
% Part 1: element-element interaction 
II = (1:NumEms);  JJ = II;
GlbMat = GlbMat + sparse(II,JJ,ArrayEE,DOFs,DOFs);
% Part 2: element-edge and edge-element interaction 
II = (1:NumEms);
for m=1:maxn
  JJ = gpos(:,m);
  GlbMat = GlbMat + sparse(II,JJ,ArrayEG(:,m),DOFs,DOFs);
  GlbMat = GlbMat + sparse(JJ,II,ArrayEG(:,m),DOFs,DOFs);
end
% Part 3: edge-edge interaction 
for i=1:maxn
  II = gpos(:,i);
  for m=1:maxn 
    JJ = gpos(:,m);
    GlbMat = GlbMat + sparse(II,JJ,ArrayGG(:,i,m),DOFs,DOFs);
  end
end

%% Assembling GlbRHS 
% Applying the chosen Gaussian quadrature for triangles 
GlbRHS = zeros(DOFs,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1); 
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  fintgrl = 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] ;
      fintgrl = fintgrl + GAUSSQUAD.TRIG(k,4) * EqnBC.fxnf(qp) * ar;
    end    
  end
  GlbRHS(ie) = fintgrl;
end

%% Incorporating boundary conditions: Neumann as natural 
% Working on all these edges simultaneously 
if NumNeumannEgs>0
  x1 = PlygnMesh.node(PlygnMesh.edge(NeumannEdge,1),1);
  y1 = PlygnMesh.node(PlygnMesh.edge(NeumannEdge,1),2);
  x2 = PlygnMesh.node(PlygnMesh.edge(NeumannEdge,2),1);
  y2 = PlygnMesh.node(PlygnMesh.edge(NeumannEdge,2),2);
  LenNeumannEg = sqrt((x2-x1).^2+(y2-y1).^2);
  % BndryLabel = PlygnMesh.BndryEdge(NeumannEdge);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    GlbRHS(NumEms+NeumannEdge) = GlbRHS(NumEms+NeumannEdge) ...
      - GAUSSQUAD.LINE(k,3) * EqnBC.fxnuN(qp) .* LenNeumannEg;
    % NEW HERE !
    % GlbRHS(NumEms+NeumannEdge) = GlbRHS(NumEms+NeumannEdge)...
    %   - GAUSSQUAD.LINE(k,3)*EqnBC.fxnuN(qp,BndryLabel).*LenNeumannEg;
  end
end

%% Handling boundary conditions: Dirichlet as essential 
% Working on all these edges simultaneously 
if NumDirichletEgs>0  
  pDavg = zeros(NumDirichletEgs,1);
  x1 = PlygnMesh.node(PlygnMesh.edge(DirichletEdge,1),1);
  y1 = PlygnMesh.node(PlygnMesh.edge(DirichletEdge,1),2);
  x2 = PlygnMesh.node(PlygnMesh.edge(DirichletEdge,2),1);
  y2 = PlygnMesh.node(PlygnMesh.edge(DirichletEdge,2),2);
  % BndryLabel = PlygnMesh.BndryEdge(DirichletEdge);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    pDavg = pDavg + GAUSSQUAD.LINE(k,3)*EqnBC.fxnpD(qp);
    % NEW HERE !
    % pDavg = pDavg + GAUSSQUAD.LINE(k,3)*EqnBC.fxnpD(qp,BndryLabel);
  end
end

%% Reducing ...
flag = zeros(DOFs,1);
flag(NumEms+DirichletEdge) = ones(NumDirichletEgs,1);
EmFreeEg = find(~flag);
sln = zeros(DOFs,1);
sln(NumEms+DirichletEdge) = pDavg;
GlbRHS = GlbRHS - GlbMat*sln;

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

return;