function [sln,GlbMat,GlbRHS] = Darcy_MFEM_RectRT0Q0s_AsmSlv(...
  EqnBC,BndryCondType,RectMesh,PermKs,GAUSSQUAD)
%% Darcy: MFEM(RT[0],Q0) on a rectangular mesh: Assembly & Solving 
% Cf. Lin,Liu,Sadre, J.Comput.Appl.Math., 2015 
% Assuming PermKs is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2018/06 

%% Mesh info 
NumEms = RectMesh.NumEms;
NumEgs = RectMesh.NumEgs;
x1 = RectMesh.node(RectMesh.elem(:,1),1);
y1 = RectMesh.node(RectMesh.elem(:,1),2);
x2 = RectMesh.node(RectMesh.elem(:,3),1);
y2 = RectMesh.node(RectMesh.elem(:,3),2);
area = (x2-x1).*(y2-y1);
AllEg = RectMesh.node(RectMesh.edge(:,2),:)...
      - RectMesh.node(RectMesh.edge(:,1),:);
LenEg = sqrt(AllEg(:,1).^2+AllEg(:,2).^2);

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

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

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

%% Darcy_MFEM_RectRT0Q0s: Assembling the global coefficient matrix 
% disp('Assembling the global coefficient matrix...'); 
% Part 1: Interaction among edge-based vec.bas.fxns. 
% For vertical edges 
II = [RectMesh.elem2edge(:,2),RectMesh.elem2edge(:,4),...
      RectMesh.elem2edge(:,2),RectMesh.elem2edge(:,4)];
JJ = [RectMesh.elem2edge(:,2),RectMesh.elem2edge(:,2),...
      RectMesh.elem2edge(:,4),RectMesh.elem2edge(:,4)];
GM = [(area/3)./PermKs, (area/6)./PermKs, (area/6)./PermKs, (area./3)./PermKs];
GlbMat = GlbMat + sparse(II(:),JJ(:),GM(:),DOFs,DOFs);
% For horizontal edges 
II = [RectMesh.elem2edge(:,1),RectMesh.elem2edge(:,3),...
      RectMesh.elem2edge(:,1),RectMesh.elem2edge(:,3)];
JJ = [RectMesh.elem2edge(:,1),RectMesh.elem2edge(:,1),...
      RectMesh.elem2edge(:,3),RectMesh.elem2edge(:,3)];
GM = [(area/3)./PermKs, (area/6)./PermKs, (area/6)./PermKs, (area./3)./PermKs];
GlbMat = GlbMat + sparse(II(:),JJ(:),GM(:),DOFs,DOFs);
% Part 2&3: Interaction among edge-based vec.bas.fxns. & elt.-based sca.bas.fxns. 
II = NumEgs + repmat((1:NumEms)',1,4);
JJ = RectMesh.elem2edge;
sn = repmat([-1,1,1,-1],NumEms,1);
pdiv = [LenEg(RectMesh.elem2edge(:,1)),LenEg(RectMesh.elem2edge(:,2)),...
        LenEg(RectMesh.elem2edge(:,3)),LenEg(RectMesh.elem2edge(:,4))].*sn;
MatB = sparse(II(:),JJ(:),pdiv(:),DOFs,DOFs);
GlbMat = GlbMat + MatB + MatB';

%% Assembling the global right-hand side: GlbRHS 
% GlbRHS = zeros(DOFs,1);
% Applying a Gaussian quadrature on rectangular elements 
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEms,2);
for k=1:NumQuadPts
  qp(:,1) = GAUSSQUAD.RECT(k,1)*x1 + GAUSSQUADRECT1(k,1)*x2;
  qp(:,2) = GAUSSQUAD.RECT(k,2)*y1 + GAUSSQUADRECT1(k,2)*y2;
  GlbRHS(NumEgs+(1:NumEms)') = GlbRHS(NumEgs+(1:NumEms)')...
    - GAUSSQUAD.RECT(k,3)*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 = RectMesh.node(RectMesh.edge(DirichletEdge,1),1);
  y1 = RectMesh.node(RectMesh.edge(DirichletEdge,1),2);
  x2 = RectMesh.node(RectMesh.edge(DirichletEdge,2),1);
  y2 = RectMesh.node(RectMesh.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

%% Adjusting the global linear system: Enforcing Neumann condition 
% disp('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 = RectMesh.node(RectMesh.edge(NeumannEdge,1),1);
  x2 = RectMesh.node(RectMesh.edge(NeumannEdge,2),1);
  y1 = RectMesh.node(RectMesh.edge(NeumannEdge,1),2);
  y2 = RectMesh.node(RectMesh.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
% Reducing...
GlbRHS = GlbRHS - GlbMat*sln;

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

return;