function [sln,MatA,MatB,GlbRHS] = Darcy_MFEM_TriRT0P0s_AsmSlv(...
  EqnBC,TriMesh,PermKs,DirichletEdge,NeumannEdge,GAUSSQUAD) 
%% Darcy: MFEM(RT0,P0) on a triangular mesh: Assembly & Solving 
% Assuming PermKs is an elementwise constant scalar 
% JL20130816: TO BE REVISED FOR EFFICIENCY 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;
% Assuming TriMesh.flag=3 
area = TriMesh.area;
LenEg = TriMesh.LenEg;
sn = TriMesh.SignEmEg;

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

%% Inverse permeability
PermKs1 = 1./PermKs;

%% 
GMRT0 = Hdiv_TriRT0_EgBas_GramMat(TriMesh);

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

%% Assembling the global coefficient matrix 
% disp('Assembling the global coefficient matrix...'); 
MatA = sparse(NumEgs,NumEgs); 
for i=1:3
  II = TriMesh.elem2edge(:,i);
  for j=1:3
    JJ = TriMesh.elem2edge(:,j);
    MatA = MatA + sparse(II,JJ,PermKs1.*GMRT0(:,i,j).*sn(:,i).*sn(:,j),...
      NumEgs,NumEgs);
  end
end
% Matrix for elementwise pressure-divergence integrals 
MatB = sparse(NumEms,NumEgs);
II = (1:NumEms)';
for j=1:3
  JJ = TriMesh.elem2edge(:,j);
  MatB = MatB + sparse(II,JJ,-LenEg(TriMesh.elem2edge(:,j)).*sn(:,j),...
    NumEms,NumEgs);
end
% Final assembly 
MatC = sparse(NumEms,NumEms);
GlbMat = [MatA,MatB'; MatB,MatC];

%% Assembling the global right-hand side: GlbRHS 
GlbRHS = zeros(DOFs,1);
% Applying a Gaussian quadrature on triangular elements 
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),:);
  % Notice the negative sign!
  GlbRHS(NumEgs+(1:NumEms)') = GlbRHS(NumEgs+(1:NumEms)')...
    - GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp).*area;
end

%% Incorporating boundary conditions: Dirichlet as natural 
% Note: only the very edge bas. fxn. assoc. w/ a Dirichlet edge is used 
% Approach I: working on all these edges simultaneously 
if NumDirichletEgs>0
  x1 = TriMesh.node(TriMesh.edge(DirichletEdge,1),1);
  y1 = TriMesh.node(TriMesh.edge(DirichletEdge,1),2);
  x2 = TriMesh.node(TriMesh.edge(DirichletEdge,2),1);
  y2 = TriMesh.node(TriMesh.edge(DirichletEdge,2),2);
  LenDirichletEg = sqrt((x2-x1).^2+(y2-y1).^2);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  qp = zeros(NumDirichletEgs,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(DirichletEdge) = GlbRHS(DirichletEdge)...
      - GAUSSQUAD.LINE(k,3)*EqnBC.fxnpD(qp).*LenDirichletEg;
    % Notice the above negative sign! 
  end
end
% % Approach II: (the old-style) loop thru these edges 
% % using the 1-pt Gaussian quadrature 
% if NumDirichletEgs>0
%   for ig=1:NumDirichletEgs
%     LblEg = DirichletEdge(ig);
%     LblVrtx = TriMesh.edge(LblEg,1:2);
%     vrtx = TriMesh.node(LblVrtx,1:2);
%     LenEg = norm(vrtx(2,:)-vrtx(1,:),2);
%     cntr = 0.5*(vrtx(1,:)+vrtx(2,:));
%     GlbRHS(LblEg) = -EqnBC.fxnpD(cntr)*LenEg;
%   end
% end

%% Adjusting the global linear system 
% For reducing... 
flag = zeros(DOFs,1);
flag(NeumannEdge) = ones(NumNeumannEgs,1);
EmFreeEg = find(~flag);
sln = zeros(DOFs,1);
% Incorporating boundary conditions: Neumann as essential 
% Approach I: working on all these edges simultaneously 
if NumNeumannEgs>0
  uNavg = zeros(NumNeumannEgs,1);
  x1 = TriMesh.node(TriMesh.edge(NeumannEdge,1),1);
  x2 = TriMesh.node(TriMesh.edge(NeumannEdge,2),1);
  y1 = TriMesh.node(TriMesh.edge(NeumannEdge,1),2);
  y2 = TriMesh.node(TriMesh.edge(NeumannEdge,2),2);
  qp = zeros(NumNeumannEgs,2);
  NumQuadPts = size(GAUSSQUAD.LINE,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;
    uNavg = uNavg + GAUSSQUAD.LINE(k,3)*EqnBC.fxnuN(qp);
  end
  sln(NeumannEdge) = uNavg;
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);
%     cntr = 0.5*(vrtx(1,:)+vrtx(2,:));
%     uNavg = EqnBC.fxnuN(cntr);  % *LenEg/LenEg
%     sln(LblEg) = uNavg;
%   end
% end
% Reducing...
GlbRHS = GlbRHS - GlbMat*sln;

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

%% JL20141218: Double-check code for reversing sign 
sln = - sln;

return;