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

%% Mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;
% Assuming TriMesh.flag>=2 
e2g = TriMesh.elem2edge;
area = TriMesh.area;

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

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

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

%% Auxiliary arrays related to discrete weak gradients across all elements 
GMK = Hdiv_TriRT0_NmlzBas_GramMatK(TriMesh,PermK);
CDWGB = WG_TriP0P0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(TriMesh);
AryEE = zeros(NumEms,1);
AryEG = zeros(NumEms,3);
AryGG = zeros(NumEms,3,3);
%% JL20151225: TO BE REVISED FOR EFFECIENCY !!
for ie=1:NumEms
  EltGM = squeeze(GMK(ie,:,:));
  AryEE(ie) = squeeze(CDWGB(ie,1,:))' * EltGM * squeeze(CDWGB(ie,1,:));
  for i=1:3
    AryEG(ie,i) = squeeze(CDWGB(ie,1,:))' * EltGM * squeeze(CDWGB(ie,1+i,:));
    for j=1:3
      AryGG(ie,i,j) = squeeze(CDWGB(ie,1+i,:))' * EltGM * squeeze(CDWGB(ie,1+j,:));
    end
  end
end

%% Assembling the global coefficient matrix GlbMat 
% In recognition of the structure of the element stiffness matrice, 
% the assembly of the global stiffness matrix is divided into 3 parts.
% Part 1: element-element interaction is diagonal b/c const in (P0,P0;RT0)
II = (1:NumEms);  JJ = II;
GlbMat = sparse(II,JJ,AryEE,DOFs,DOFs);
% Part 2: element-edge interaction 
II = (1:NumEms);
for j=1:3 
  JJ = NumEms + e2g(:,j);
  GlbMat = GlbMat + sparse(II,JJ,AryEG(:,j),DOFs,DOFs);
  GlbMat = GlbMat + sparse(JJ,II,AryEG(:,j),DOFs,DOFs);
end
% Part 3: edge-edge interaction 
for i=1:3
  II = NumEms + e2g(:,i);
  for j=i:3  % Utilizing symmetry 
    JJ = NumEms + e2g(:,j);
    if (j==i) 
      GlbMat = GlbMat + sparse(II,JJ,AryGG(:,i,i),DOFs,DOFs);
    else
      GlbMat = GlbMat + sparse(II,JJ,AryGG(:,i,j),DOFs,DOFs);
      GlbMat = GlbMat + sparse(JJ,II,AryGG(:,i,j),DOFs,DOFs);
    end
  end
end

%% Assembling GlbRHS 
% Applying a Gaussian quadrature on triangular elements 
GlbRHS = zeros(DOFs,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * TriMesh.node(TriMesh.elem(:,1),:)...
     + GAUSSQUAD.TRIG(k,2) * TriMesh.node(TriMesh.elem(:,2),:)...
     + GAUSSQUAD.TRIG(k,3) * TriMesh.node(TriMesh.elem(:,3),:);
  GlbRHS(1:NumEms) = GlbRHS(1:NumEms) + GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp);
end
GlbRHS(1:NumEms) = GlbRHS(1:NumEms).*area;
% % The old-style code 
% for ie=1:NumEms
%   LblVrtx = TriMesh.elem(ie,1:3);
%   vrtx = TriMesh.node(LblVrtx,1:2);
%   QuadPt = GAUSSQUAD.TRIG(:,1:3)*vrtx;
%   fval = EqnBC.fxnf(QuadPt);
%   GlbRHS(ie) = dot(fval,GAUSSQUAD.TRIG(:,4))*area(ie);  
% end

%% Incorporating boundary conditions: Neumann as natural 
% Approach I: working on all these edges simultaneously 
if NumNeumannEgs>0
  x1 = TriMesh.node(TriMesh.edge(NeumannEdge,1),1);
  y1 = TriMesh.node(TriMesh.edge(NeumannEdge,1),2);
  x2 = TriMesh.node(TriMesh.edge(NeumannEdge,2),1);
  y2 = TriMesh.node(TriMesh.edge(NeumannEdge,2),2);
  LenNeumannEg = sqrt((x2-x1).^2+(y2-y1).^2);
  % LenNeumannEg = norm(TriMesh.node(TriMesh.edge(NeumannEdge,2),:)...
  %   - TriMesh.node(TriMesh.edge(NeumannEdge,1),:),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
% % Approach II: (the old-style) loop thru these edges 
% if NumNeumannEgs>0 
%   for ig=1:NumNeumannEgs
%     LblEg = NeumannEdge(ig);
%     LblVrtx = TriMesh.edge(LblEg,1:2);
%     vrtx = TriMesh.node(LblVrtx,1:2);
%     LenEg = norm(vrtx(2,:)-vrtx(1,:),2);
%     QuadPt = GAUSSQUAD.LINE(:,1:2)*vrtx;
%     uNval = EqnBC.fxnuN(QuadPt);
%     pos = NumEms + LblEg;
%     GlbRHS(pos) = GlbRHS(pos) - dot(uNval,GAUSSQUAD.LINE(:,3))*LenEg;
%   end
% end

%% For reducing... 
% disp('Adjusting the global linear system...'); 
flag = zeros(DOFs,1);
flag(NumEms+DirichletEdge) = ones(NumDirichletEgs,1);
EmFreeEg = find(~flag);
% sln vector to be used 
sln = zeros(DOFs,1);

%% Incorporating boundary conditions: Dirichlet as essential 
% Approach I: working on all these edges simultaneously 
% Assuming NumDirichletEgs>0
x1 = TriMesh.node(TriMesh.edge(DirichletEdge,1),1);
x2 = TriMesh.node(TriMesh.edge(DirichletEdge,2),1);
y1 = TriMesh.node(TriMesh.edge(DirichletEdge,1),2);
y2 = TriMesh.node(TriMesh.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(NumEms+DirichletEdge) = pDavg;
% % Approach II: (the old-style) loop thru these edges 
% Assuming NumDirichletEgs>0
% for ig=1:NumDirichletEgs
%   LblEg = DirichletEdge(ig);
%   LblVrtx = TriMesh.edge(LblEg,1:2);
%   vrtx = TriMesh.node(LblVrtx,1:2);
%   QuadPt = GAUSSQUAD.LINE(:,1:2)*vrtx;
%   pDval = EqnBC.fxnpD(QuadPt);
%   pDavg = dot(GAUSSQUAD.LINE(:,3),pDval);  % *LenEg/LenEg
%   sln(NumEms+LblEg) = pDavg;
% end

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

%% JL20130725: TO BE REVISED USING SCHUR COMPLEMENT 
%% Solving the reduced global linear system directly 
% disp('Solving the reduced global linear system directly...'); 
sln(EmFreeEg) = GlbMat(EmFreeEg,EmFreeEg) \ GlbRHS(EmFreeEg);

return;