function [sln,MatA,MatB,GlbRHS] = Darcy_MFEM_TriBDM1P0s_AsmSlv(...
  EqnBC, BndryCondType, TriMesh, PermKs, GAUSSQUAD) 
%% Darcy: MFEM(BDM1,P0) on a triangular mesh 
% 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=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 
NumDiriEgs = size(DirichletEdge,1);
NumNeumEgs = size(NeumannEdge,1);

%% Auxiliary quantities 
k1 = TriMesh.elem(:,1);  k2 = TriMesh.elem(:,2);  k3 = TriMesh.elem(:,3);
x1 = TriMesh.node(k1,1);  y1 = TriMesh.node(k1,2); 
x2 = TriMesh.node(k2,1);  y2 = TriMesh.node(k2,2);
x3 = TriMesh.node(k3,1);  y3 = TriMesh.node(k3,2);
eg1 = [x3,y3] - [x2,y2];
eg2 = [x1,y1] - [x3,y3];
eg3 = [x2,y2] - [x1,y1];
lg1 = sqrt(eg1(:,1).^2+eg1(:,2).^2);
lg2 = sqrt(eg2(:,1).^2+eg2(:,2).^2);
lg3 = sqrt(eg3(:,1).^2+eg3(:,2).^2);
tan = zeros(NumEms,3,2);
tan(:,1,:) = eg1./[lg1,lg1];
tan(:,2,:) = eg2./[lg2,lg2];
tan(:,3,:) = eg3./[lg3,lg3];
% nml = zeros(NumEms,3,2);
% nml(:,1,:) = [tan(:,1,2),-tan(:,1,1)];
% nml(:,2,:) = [tan(:,2,2),-tan(:,2,1)];
% nml(:,3,:) = [tan(:,3,2),-tan(:,3,1)];
% cosine theta 
ct1 = (lg2.^2+lg3.^2-lg1.^2)./(2*lg2.*lg3);
ct2 = (lg3.^2+lg1.^2-lg2.^2)./(2*lg3.*lg1);
ct3 = (lg1.^2+lg2.^2-lg3.^2)./(2*lg1.*lg2);
% sine theta 
st1 = sqrt(1-ct1.^2);
st2 = sqrt(1-ct2.^2);
st3 = sqrt(1-ct3.^2);
% Dot products among the three edge tangetial vectors for all elements 
dptt = zeros(NumEms,3,3);
for i=1:3
  for j=1:3
    dptt(:,i,j) = dot(tan(:,i,:),tan(:,j,:),3);
  end
end

%% Inverse permeability
PermKs1 = 1./PermKs;

%% 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);

%% Gram matrices in BDM1 eg.bas.bas.fxns.(local) 
% For each element, the Gram matrix is 6x6 
GMBDM1 = zeros(NumEms,6,6);
% diagonal entries 
GMBDM1(:,1,1) = (area/6)./(st2.^2);  GMBDM1(:,2,2) = (area/6)./(st3.^2);
GMBDM1(:,3,3) = (area/6)./(st3.^2);  GMBDM1(:,4,4) = (area/6)./(st1.^2);
GMBDM1(:,5,5) = (area/6)./(st1.^2);  GMBDM1(:,6,6) = (area/6)./(st2.^2);
% 1st row off-diagonals 
GMBDM1(:,1,2) = -dptt(:,3,2)./(st2.*st3).*(area/12);
GMBDM1(:,1,3) =  dptt(:,3,1)./(st2.*st3).*(area/12);
GMBDM1(:,1,4) = -dptt(:,3,3)./(st2.*st1).*(area/12);
GMBDM1(:,1,5) =  dptt(:,3,2)./(st2.*st1).*(area/12);
GMBDM1(:,1,6) = -dptt(:,3,1)./(st2.*st2).*(area/6);
% 2nd row off-diagonals 
GMBDM1(:,2,3) = -dptt(:,2,1)./(st3.*st3).*(area/6);
GMBDM1(:,2,4) =  dptt(:,2,3)./(st3.*st1).*(area/12);
GMBDM1(:,2,5) = -dptt(:,2,2)./(st3.*st1).*(area/12);
GMBDM1(:,2,6) =  dptt(:,2,1)./(st3.*st2).*(area/12);
% 3rd row off-diagonals 
GMBDM1(:,3,4) = -dptt(:,1,3)./(st3.*st1).*(area/12);
GMBDM1(:,3,5) =  dptt(:,1,2)./(st3.*st1).*(area/12);
GMBDM1(:,3,6) = -dptt(:,1,1)./(st3.*st2).*(area/12);
% 4th row off-diagonals 
GMBDM1(:,4,5) = -dptt(:,3,2)./(st1.*st1).*(area/6);
GMBDM1(:,4,6) =  dptt(:,3,1)./(st1.*st2).*(area/12);
% 5th row off-diagonals 
GMBDM1(:,5,6) = -dptt(:,2,1)./(st1.*st2).*(area/12);
% symmetry 
for i=2:6 
  for j=1:(i-1) 
    GMBDM1(:,i,j) = GMBDM1(:,j,i);
  end
end
% Adjusting 
for ie=1:NumEms
  for j=1:3
    if (TriMesh.SignEmEg(ie,j)==-1) 
      GMBDM1(ie,:,:) = MatPermu{j}*squeeze(GMBDM1(ie,:,:))*MatPermu{j};
    end
  end
end

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

%% Assembling the global coefficient matrix 
% disp('Assembling the global coefficient matrix...'); 
MatA = sparse(2*NumEgs,2*NumEgs); 
% JL20131001: TO BE REVISED FOR EFFICIENCY 
for i=1:3
  II = 2*TriMesh.elem2edge(:,i)-1;
  for j=1:3
    JJ = 2*TriMesh.elem2edge(:,j)-1;
    MatA = MatA + sparse(II,JJ,PermKs1.*GMBDM1(:,2*i-1,2*j-1).*sn(:,i).*sn(:,j),...
      2*NumEgs,2*NumEgs);
    JJ = 2*TriMesh.elem2edge(:,j);  
    MatA = MatA + sparse(II,JJ,PermKs1.*GMBDM1(:,2*i-1,2*j).*sn(:,i).*sn(:,j),...
      2*NumEgs,2*NumEgs);
  end
  II = 2*TriMesh.elem2edge(:,i);
  for j=1:3
    JJ = 2*TriMesh.elem2edge(:,j)-1;
    MatA = MatA + sparse(II,JJ,PermKs1.*GMBDM1(:,2*i,2*j-1).*sn(:,i).*sn(:,j),...
      2*NumEgs,2*NumEgs);
    JJ = 2*TriMesh.elem2edge(:,j);  
    MatA = MatA + sparse(II,JJ,PermKs1.*GMBDM1(:,2*i,2*j).*sn(:,i).*sn(:,j),...
      2*NumEgs,2*NumEgs);
  end
end
% Gradients of the Lagrangian P1 bas.fxns. = Barycentric coordinates 
% CP1BasGrad = zeros(NumEms,3,2);  % All elts., 3 bas.fxns.
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);
% Divergence of the BDM1 bas.fxns.
DivBDM1Bas = zeros(NumEms,6);  % all elts., 3 edges, 3 ends 
DivBDM1Bas(:,1) =  dot(CP1BasGrad(:,2,:),tan(:,3,:),3)./st2;
DivBDM1Bas(:,2) = -dot(CP1BasGrad(:,3,:),tan(:,2,:),3)./st3;
DivBDM1Bas(:,3) =  dot(CP1BasGrad(:,3,:),tan(:,1,:),3)./st3;
DivBDM1Bas(:,4) = -dot(CP1BasGrad(:,1,:),tan(:,3,:),3)./st1;
DivBDM1Bas(:,5) =  dot(CP1BasGrad(:,1,:),tan(:,2,:),3)./st1;
DivBDM1Bas(:,6) = -dot(CP1BasGrad(:,2,:),tan(:,1,:),3)./st2;
% Matrix for elementwise pressure-divergence integrals 
pdse = zeros(NumEms,6);
for j=1:6 
  pdse(:,j) = -DivBDM1Bas(:,j).*area;
end
% Adjusting 
for ie=1:NumEms
  for j=1:3
    if (TriMesh.SignEmEg(ie,j)==-1) 
      pdse(ie,1:6) = pdse(ie,:)*MatPermu{j};
    end
  end
end
% Assembling MatB 
MatB = sparse(NumEms,2*NumEgs);
II = (1:NumEms);
for j=1:3
  JJ = 2*TriMesh.elem2edge(:,j)-1;
  MatB = MatB + sparse(II,JJ,pdse(:,2*j-1).*sn(:,j),NumEms,2*NumEgs);
  JJ = 2*TriMesh.elem2edge(:,j);
  MatB = MatB + sparse(II,JJ,pdse(:,2*j).*sn(:,j),NumEms,2*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 
% qp = zeros(NumEms,2);
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(2*NumEgs+(1:NumEms)) = GlbRHS(2*NumEgs+(1:NumEms))...
    - GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp).*area;  % Notice "-" 
end
% % JL20131005: Old-style 
% GlbRHS = zeros(DOFs,1);
% NumQuadPts = size(GAUSSQUAD.TRIG,1);
% for ie=1:NumEms
%   EltRHS = 0;
%   for k=1:NumQuadPts
%     qp = GAUSSQUAD.TRIG(k,1) * TriMesh.node(TriMesh.elem(ie,1),:)...
%        + GAUSSQUAD.TRIG(k,2) * TriMesh.node(TriMesh.elem(ie,2),:)...
%        + GAUSSQUAD.TRIG(k,3) * TriMesh.node(TriMesh.elem(ie,3),:);
%     EltRHS = EltRHS + GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp)*area(ie);
%   end
%   GlbRHS(2*NumEgs+ie) = - EltRHS;  % Notice "-" 
% end

%% Incorporating boundary conditions: Dirichlet as natural 
% Note: only the very edge bas. fxn. assoc. w/ a Dirichlet edge is used 
% Approach II: Old-style: running thru all edges 
DiriRHS = zeros(NumDiriEgs,2);  % Assuming NumDiriEgs>0 
NumQuadPts = size(GAUSSQUAD.LINE,1);
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);
  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(ig,1) = DiriRHS(ig,1) + GAUSSQUAD.LINE(k,3)*(pD.*t1);
    DiriRHS(ig,2) = DiriRHS(ig,2) + GAUSSQUAD.LINE(k,3)*(pD.*t2);
  end
  DiriRHS(ig,1) = DiriRHS(ig,1).*LenDiriEg;
  DiriRHS(ig,2) = DiriRHS(ig,2).*LenDiriEg;
  % Incorporating: noting the negative signs below! 
  GlbRHS(2*lbleg-1) = GlbRHS(2*lbleg-1) - DiriRHS(ig,1);
  GlbRHS(2*lbleg) = GlbRHS(2*lbleg) - DiriRHS(ig,2);
end

%% Adjusting the global linear system 
% disp('Adjusting the global linear system...'); 
% For reducing... 
flag = ones(DOFs,1);
flag(2*NeumannEdge-1) = zeros(NumNeumEgs,1);
flag(2*NeumannEdge) = zeros(NumNeumEgs,1);
EmFreeEg = find(flag);
sln = zeros(DOFs,1);
% Incorporating boundary conditions: Neumann as essential 
% 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;
    % Solving 
    uNproj = [4,-2;-2,4]*NeumRHS;
    % Enforcing 
    sln(2*lbleg-1) = uNproj(1);
    sln(2*lbleg) = uNproj(2);
  end
end
% Reducing...
GlbRHS = GlbRHS - GlbMat*sln;

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

return;