function [sln] = Darcy_MFEM_TriRT0P0_AsmSlv(...
  EqnBC,BndryCondType,TriMesh,PermK,GAUSSQUAD) 
%% Darcy: MFEM(RT0,P0): Assembly & Solving for a triangular mesh 
% Assuming PermK is an elementwise constant 2*2 SPD matrix 
% James Liu, ColoState; 2012/07--2021/01 

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

%% Sorting out boundary edges: Dirichlet, 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 = NumEgs + NumEms;
% GlbMat = sparse(DOFs,DOFs);
% GlbRHS = zeros(DOFs,1); 

%% Elementwise Gram matrix edge-based basis functions with weight K^{-1} 
GMRT0K1 = Hdiv_TriRT0_EgBas_GramMatK1(TriMesh,PermK);

%% Assembling the global "stiffness" matrix 
% MatA for global Gram matrix edge-based basis functions with weight K^{-1}
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, GMRT0K1(:,i,j).*sn(:,i).*sn(:,j),...
      NumEgs,NumEgs);
  end
end
% MatB for global matrix for pressure-divergence interaction 
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 
% MatC is 0/empty for now 
MatC = sparse(NumEms,NumEms); 
% Final assembly 
GlbMat = [MatA,MatB'; MatB,MatC];

%% Assembling the global right-hand side: GlbRHS 
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(NumEgs+(1:NumEms)') = GlbRHS(NumEgs+(1:NumEms)')...
    - GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp).*area;  % NOTE: Negative sign 
end

%% Incorporating boundary conditions: Dirichlet as natural 
% Note: only the very edge bas. fxn. assoc. w/ a Dirichlet edge is used 
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);
  qp = zeros(NumDirichletEgs,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;
    GlbRHS(DirichletEdge) = GlbRHS(DirichletEdge)...
      - GAUSSQUAD.LINE(k,3)*EqnBC.fxnpD(qp).*LenDirichletEg;  % NOTE: "-" 
  end
end

%% For reducing 
flag = zeros(DOFs,1);
flag(NeumannEdge) = ones(NumNeumannEgs,1);
sln = zeros(DOFs,1);

%% Incorporating boundary conditions: Neumann as essential 
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

%% Reducing 
GlbRHS = GlbRHS - GlbMat*sln;

%% Solving the reduced global linear system directly 
% disp('Solving the reduced global linear system directly...'); 
EmFreeEg = find(~flag);
sln(EmFreeEg) = GlbMat(EmFreeEg,EmFreeEg)\GlbRHS(EmFreeEg);

return;