function [sln,GlbMat,GlbRHS] = Darcy_WG_TriP0P1BDM1s_AsmSlv(...
  EqnBC,BndryCondType,TriMesh,PermKs,GAUSSQUAD)
%% Darcy: WG(P0,P1;BDM1) on a triangular mesh: Assembly & Solving 
% Assuming PermKs is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2018/06 

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

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

%% Boundary info 
NumDiriEgs = size(DirichletEdge,1);
NumNeumEgs = size(NeumannEdge,1);

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

%% Computing coeff. disc.wk.grad. in BDM1 eg.bas.fxns. of WG bas.fxns.
CDWGB = WG_TriP0P1BDM1_CofBDM1EgBas_DiscWkGradBasFxn(TriMesh);

%% Computing Gram matrices of BDM1 edge-based bas.fxns. on a triangular mesh 
GMBDM1 = Hdiv_TriBDM1_EgBas_GramMat(TriMesh);

%% Assembling the global coefficient matrix GlbMat 
% In recognition of the structure of the element stiffness matrice, 
% entries & assembly of GlbMat are divided into 3 parts.
EntryE = zeros(NumEms,1);  
EntryEG = zeros(NumEms,6);
EntryGG = zeros(NumEms,6,6);
for ie=1:NumEms
  EltGM = squeeze(GMBDM1(ie,1:6,1:6));
  EntryE(ie) = CDWGB(ie,1:6,7)*(EltGM*CDWGB(ie,1:6,7)');
  for i=1:6
    EntryEG(ie,i) = CDWGB(ie,1:6,7)*(EltGM*CDWGB(ie,1:6,i)');
    for j=1:6
      EntryGG(ie,i,j) = CDWGB(ie,1:6,i)*(EltGM*CDWGB(ie,1:6,j)');
    end
  end
end

%% Adjusting entries
% Permutation matrices for adjusting 
MatPermu = cell(3,1);
for i=1:3; MatPermu{i} = zeros(6);  end
MatPermu{1}(1,2)=1;  MatPermu{1}(2,1)=1;  
MatPermu{1}(3:4,3:4)=eye(2);  MatPermu{1}(5:6,5:6)=eye(2);
MatPermu{2}(3,4)=1;  MatPermu{2}(4,3)=1;  
MatPermu{2}(1:2,1:2)=eye(2);  MatPermu{2}(5:6,5:6)=eye(2);
MatPermu{3}(5,6)=1;  MatPermu{3}(6,5)=1;
MatPermu{3}(1:2,1:2)=eye(2);  MatPermu{3}(3:4,3:4)=eye(2);
for ie=1:NumEms
  for j=1:3
    if TriMesh.SignEmEg(ie,j)==-1 
      EntryEG(ie,:) = EntryEG(ie,:)*MatPermu{j};
      EntryGG(ie,:,:) = MatPermu{j}*squeeze(EntryGG(ie,:,:))*MatPermu{j};
    end
  end
end

% Part I: element-element interaction: diagonal b/c const in (P0,P1,BDM1)
II = (1:NumEms);  JJ = II; 
GlbMat = sparse(II,JJ,EntryE.*PermKs,DOFs,DOFs); 
% Part II: element-edge interaction 
II = (1:NumEms);
for j=1:3 
  JJ = NumEms + 2*TriMesh.elem2edge(:,j)-1;
  GlbMat = GlbMat + sparse(II,JJ,EntryEG(:,2*j-1).*PermKs,DOFs,DOFs);
  GlbMat = GlbMat + sparse(JJ,II,EntryEG(:,2*j-1).*PermKs,DOFs,DOFs);
  JJ = NumEms + 2*TriMesh.elem2edge(:,j);
  GlbMat = GlbMat + sparse(II,JJ,EntryEG(:,2*j).*PermKs,DOFs,DOFs);
  GlbMat = GlbMat + sparse(JJ,II,EntryEG(:,2*j).*PermKs,DOFs,DOFs);
end
% Part III: edge-edge interaction 
for i=1:3
  II = NumEms + 2*TriMesh.elem2edge(:,i)-1;
  for j=1:3
    JJ = NumEms + 2*TriMesh.elem2edge(:,j)-1;
    GlbMat = GlbMat + sparse(II,JJ,EntryGG(:,2*i-1,2*j-1).*PermKs,...
      DOFs,DOFs);
    JJ = NumEms + 2*TriMesh.elem2edge(:,j);
    GlbMat = GlbMat + sparse(II,JJ,EntryGG(:,2*i-1,2*j).*PermKs,...
      DOFs,DOFs);
  end
  II = NumEms + 2*TriMesh.elem2edge(:,i);
  for j=1:3
    JJ = NumEms + 2*TriMesh.elem2edge(:,j)-1;
    GlbMat = GlbMat + sparse(II,JJ,EntryGG(:,2*i,2*j-1).*PermKs,...
      DOFs,DOFs);
    JJ = NumEms + 2*TriMesh.elem2edge(:,j);
    GlbMat = GlbMat + sparse(II,JJ,EntryGG(:,2*i,2*j).*PermKs,...
      DOFs,DOFs);
  end
end

%% Assembling 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),:);
  GlbRHS(1:NumEms) = GlbRHS(1:NumEms)...
    + GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp).*area;
end

%% Incorporating boundary conditions: Neumann as natural 
% Approach II: Old-style: running thru all edges 
NumQuadPts = size(GAUSSQUAD.LINE,1);
if NumNeumEgs>0
  for ig=1:NumNeumEgs 
    % Edge info 
    lbleg = NeumannEdge(ig);
    ie1 = TriMesh.edge2elem(lbleg,1);
    j1 = TriMesh.WhichEdge(lbleg,1);
    js=j1+1;  if (js>3) js=mod(js,3); end 
    je=j1+2;  if (je>3) je=mod(je,3); end
    x1 = TriMesh.node(TriMesh.elem(ie1,js),1);
    y1 = TriMesh.node(TriMesh.elem(ie1,js),2);
    x2 = TriMesh.node(TriMesh.elem(ie1,je),1);
    y2 = TriMesh.node(TriMesh.elem(ie1,je),2);
    LenNeumEg = sqrt((x2-x1).^2+(y2-y1).^2);
    NeumRHS = zeros(2,1);
    for k=1:NumQuadPts
      qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
      uN = EqnBC.fxnuN(qp);
      t1 = GAUSSQUAD.LINE(k,1);
      t2 = GAUSSQUAD.LINE(k,2);
      NeumRHS(1) = NeumRHS(1) + GAUSSQUAD.LINE(k,3)*(uN*t1);
      NeumRHS(2) = NeumRHS(2) + GAUSSQUAD.LINE(k,3)*(uN*t2);
    end 
    NeumRHS = NeumRHS*LenNeumEg;
    % Incorporating: noting the negative signs below! 
    GlbRHS(NumEms+2*lbleg-1) = GlbRHS(NumEms+2*lbleg-1) - NeumRHS(1);
    GlbRHS(NumEms+2*lbleg)   = GlbRHS(NumEms+2*lbleg)   - NeumRHS(2);
  end
end

%% For reducing... 
% disp('Adjusting the global linear system...'); 
flag = ones(DOFs,1);
flag(NumEms+2*DirichletEdge-1) = zeros(NumDiriEgs,1);
flag(NumEms+2*DirichletEdge)   = zeros(NumDiriEgs,1);
EmFreeEg = find(flag);
sln = zeros(DOFs,1);
% Incorporating boundary conditions: Dirichlet as essential 
% Approach II: Old-style: running thru all edges 
% Assuming NumDiriEgs>0
NumQuadPts = size(GAUSSQUAD.LINE,1);
if NumDiriEgs>0
  for ig=1:NumDiriEgs 
    % Edge info 
    lbleg = DirichletEdge(ig);
    ie1 = TriMesh.edge2elem(lbleg,1);
    j1 = TriMesh.WhichEdge(lbleg,1);
    js=j1+1;  if (js>3) js = mod(js,3); end 
    je=j1+2;  if (je>3) je = mod(je,3); end
    x1 = TriMesh.node(TriMesh.elem(ie1,js),1);
    y1 = TriMesh.node(TriMesh.elem(ie1,js),2);
    x2 = TriMesh.node(TriMesh.elem(ie1,je),1);
    y2 = TriMesh.node(TriMesh.elem(ie1,je),2);
    % LenDiriEg = sqrt((x2-x1).^2+(y2-y1).^2);
    DiriRHS = zeros(2,1);
    for k=1:NumQuadPts
      qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
      pD = EqnBC.fxnpD(qp);
      t1 = GAUSSQUAD.LINE(k,1);
      t2 = GAUSSQUAD.LINE(k,2);
      DiriRHS(1) = DiriRHS(1) + GAUSSQUAD.LINE(k,3)*(pD*t1);
      DiriRHS(2) = DiriRHS(2) + GAUSSQUAD.LINE(k,3)*(pD*t2);
    end 
    % DiriRHS = DiriRHS*LenDiriEg;
    % Solving 
    pDproj = [4,-2;-2,4]*DiriRHS;
    % Enforcing 
    sln(NumEms+2*lbleg-1) = pDproj(1);
    sln(NumEms+2*lbleg) = pDproj(2);
  end
end
%% Reducing... 
GlbRHS = GlbRHS - GlbMat*sln;

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

return;