function [sln,GlbMat,GlbRHS] = Darcy_WG_HybridPQ0RT0_AsmSlv(...
  EqnBC,BndryCondType,HybridMesh,PermK,GAUSSQUAD)

%% Darcy: Hybrid mesh: Assembly & Solving by WG(Q0,Q0;RT[0])quadri+WG(P0,P0;RT0)
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2018/06

%% Mesh info 
NumEms  = HybridMesh.NumEms;
NumEmsQ = HybridMesh.NumEmsQ;
NumEmsT = HybridMesh.NumEmsT;
NumEgs  = HybridMesh.NumEgs;

%% Sorting out boundary edges: Dirichlet, Neumann 
DirichletEdge = find(BndryCondType(HybridMesh.BndryEdge+1)==1);
NeumannEdge   = find(BndryCondType(HybridMesh.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);

%% Assembling the global coefficient matrix GlbMat 
% GMKQ = zeros(NumEmsQ,4,4);
% GMKT = zeros(NumEmsQ,3,3);
% CDWGBQ = zeros(NumEmsQ,5,4);
% CDWGBT = zeros(NumEmsT,4,3);
[GMKQ,GMKT] = Hdiv_HybridRT0_NmlzBas_GramMatK(HybridMesh,PermK);
[CDWGBQ,CDWGBT] = WG_HybridPQ0_CofRT0NmlzBas_DiscWkGradBasFxn(HybridMesh);
% GlbMat = sparse(DOFs,DOFs);
% For the quadrilateral elements: Computing and Assembling 
ArrayEE = zeros(NumEmsQ,1);
ArrayEG = zeros(NumEmsQ,4);
ArrayGG = zeros(NumEmsQ,4,4);
%% JL20160716: TO BE REVISED FOR EFFECIENCY !!
for ie=1:NumEmsQ 
  EGM = squeeze(GMKQ(ie,:,:));
  ArrayEE(ie) = squeeze(CDWGBQ(ie,1,:))' * EGM * squeeze(CDWGBQ(ie,1,:));
  for i=1:4
    ArrayEG(ie,i) = squeeze(CDWGBQ(ie,1,:))' * EGM * squeeze(CDWGBQ(ie,1+i,:));
    for j=1:4
      ArrayGG(ie,i,j) = squeeze(CDWGBQ(ie,1+i,:))' * EGM * squeeze(CDWGBQ(ie,1+j,:));
    end
  end
end
GlbMat = sparse(DOFs,DOFs);
% Part Q.1: element-element interaction 
II = HybridMesh.labelemQ;  JJ = II;
GlbMat = GlbMat + sparse(II,JJ,ArrayEE,DOFs,DOFs);
% Part Q.2: element-edge interaction 
II = HybridMesh.labelemQ;
for j=1:4
  JJ = NumEms + HybridMesh.elemQedge(:,j);
  GlbMat = GlbMat + sparse(II,JJ,ArrayEG(:,j),DOFs,DOFs);
  GlbMat = GlbMat + sparse(JJ,II,ArrayEG(:,j),DOFs,DOFs);
end
% Part Q.3: edge-edge interaction 
for i=1:4
  II = NumEms + HybridMesh.elemQedge(:,i);
  for j=i:4  % Utilizing symmetry 
    JJ = NumEms + HybridMesh.elemQedge(:,j);
    if (j==i) 
      GlbMat = GlbMat + sparse(II,JJ,ArrayGG(:,i,j),DOFs,DOFs);
    else
      GlbMat = GlbMat + sparse(II,JJ,ArrayGG(:,i,j),DOFs,DOFs);
      GlbMat = GlbMat + sparse(JJ,II,ArrayGG(:,i,j),DOFs,DOFs);
    end
  end
end
% For the triangular elements: Computing and Assembling  
ArrayEE = zeros(NumEmsT,1);
ArrayEG = zeros(NumEmsT,3);
ArrayGG = zeros(NumEmsT,3,3);
%% JL20160716: TO BE REVISED FOR EFFECIENCY !!
for ie=1:NumEmsT 
  EGM = squeeze(GMKT(ie,:,:));
  ArrayEE(ie) = squeeze(CDWGBT(ie,1,:))' * EGM * squeeze(CDWGBT(ie,1,:));
  for i=1:3
    ArrayEG(ie,i) = squeeze(CDWGBT(ie,1,:))' * EGM * squeeze(CDWGBT(ie,1+i,:));
    for j=1:3
      ArrayGG(ie,i,j) = squeeze(CDWGBT(ie,1+i,:))' * EGM * squeeze(CDWGBT(ie,1+j,:));
    end
  end
end
% Part T.1: element-element interaction 
II = HybridMesh.labelemT;  JJ = II;
GlbMat = GlbMat + sparse(II,JJ,ArrayEE,DOFs,DOFs);
% Part T.2: element-edge interaction 
II = HybridMesh.labelemT;
for j=1:3
  JJ = NumEms + HybridMesh.elemTedge(:,j);
  GlbMat = GlbMat + sparse(II,JJ,ArrayEG(:,j),DOFs,DOFs);
  GlbMat = GlbMat + sparse(JJ,II,ArrayEG(:,j),DOFs,DOFs);
end
% Part T.3: edge-edge interaction 
for i=1:3
  II = NumEms + HybridMesh.elemTedge(:,i);
  for j=i:3  % Utilizing symmetry 
    JJ = NumEms + HybridMesh.elemTedge(:,j);
    if (j==i) 
      GlbMat = GlbMat + sparse(II,JJ,ArrayGG(:,i,j),DOFs,DOFs);
    else
      GlbMat = GlbMat + sparse(II,JJ,ArrayGG(:,i,j),DOFs,DOFs);
      GlbMat = GlbMat + sparse(JJ,II,ArrayGG(:,i,j),DOFs,DOFs);
    end
  end
end

%% Computing and Assembling GlbRHS 
GlbRHS = zeros(DOFs,1);
% For the quadrilateral elements 
LblQ = HybridMesh.labelemQ;
k1 = HybridMesh.elemQ(:,1);  k2 = HybridMesh.elemQ(:,2);  
k3 = HybridMesh.elemQ(:,3);  k4 = HybridMesh.elemQ(:,4);
x1 = HybridMesh.node(k1,1);  y1 = HybridMesh.node(k1,2);
x2 = HybridMesh.node(k2,1);  y2 = HybridMesh.node(k2,2);
x3 = HybridMesh.node(k3,1);  y3 = HybridMesh.node(k3,2);
x4 = HybridMesh.node(k4,1);  y4 = HybridMesh.node(k4,2);
CofA = zeros(NumEmsQ,4);      CofB = zeros(NumEmsQ,4);
CofA(:,1) = x1;               CofB(:,1) = y1;
CofA(:,2) = x2-x1;            CofB(:,2) = y2-y1;
CofA(:,3) = x4-x1;            CofB(:,3) = y4-y1;
CofA(:,4) = (x1+x3)-(x2+x4);  CofB(:,4) = (y1+y3)-(y2+y4);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEmsQ,2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp(:,1) = CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat;
  qp(:,2) = CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat;
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  GlbRHS(LblQ) = GlbRHS(LblQ)...
    + GAUSSQUAD.RECT(k,3) * jac .* EqnBC.fxnf(qp);
end
% For the triangular elements 
LblT = HybridMesh.labelemT;
k1 = HybridMesh.elemT(:,1);  
k2 = HybridMesh.elemT(:,2);  
k3 = HybridMesh.elemT(:,3);
x1 = HybridMesh.node(k1,1);  y1 = HybridMesh.node(k1,2); 
x2 = HybridMesh.node(k2,1);  y2 = HybridMesh.node(k2,2);
x3 = HybridMesh.node(k3,1);  y3 = HybridMesh.node(k3,2);
area = 0.5*((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * HybridMesh.node(HybridMesh.elemT(:,1),:)...
     + GAUSSQUAD.TRIG(k,2) * HybridMesh.node(HybridMesh.elemT(:,2),:)...
     + GAUSSQUAD.TRIG(k,3) * HybridMesh.node(HybridMesh.elemT(:,3),:);
  GlbRHS(LblT) = GlbRHS(LblT) + GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp);
end
GlbRHS(LblT) = GlbRHS(LblT).*area;

%% Incorporating boundary conditions: Neumann as natural 
% Approach I: working on all these edges simultaneously 
if NumNeumannEgs>0
  x1 = HybridMesh.node(HybridMesh.edge(NeumannEdge,1),1);
  y1 = HybridMesh.node(HybridMesh.edge(NeumannEdge,1),2);
  x2 = HybridMesh.node(HybridMesh.edge(NeumannEdge,2),1);
  y2 = HybridMesh.node(HybridMesh.edge(NeumannEdge,2),2);
  LenNeumannEg = sqrt((x2-x1).^2+(y2-y1).^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

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

%% Incorporating boundary conditions: Dirichlet as essential 
% Approach I: working on all these edges simultaneously 
if NumDirichletEgs>0 
  x1 = HybridMesh.node(HybridMesh.edge(DirichletEdge,1),1);
  y1 = HybridMesh.node(HybridMesh.edge(DirichletEdge,1),2);
  x2 = HybridMesh.node(HybridMesh.edge(DirichletEdge,2),1);
  y2 = HybridMesh.node(HybridMesh.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;
end

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

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

return;