function [sln,GlbMat,GlbRHS] = Darcy_FVMD_TriP1s_AsmSlv(...
  EqnBC,BndryCondType,TriMesh,PermKs,GAUSSQUAD)
%% Darcy: FVM node-oriented: Tri.mesh: P1 
% Assuming PermKs is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2018/06 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumNds = TriMesh.NumNds;

%% For the dual mesh 
% Getting all 7 vertices 
vrtx1 = TriMesh.node(TriMesh.elem(:,1),:);
vrtx2 = TriMesh.node(TriMesh.elem(:,2),:);
vrtx3 = TriMesh.node(TriMesh.elem(:,3),:);
vrtx4 = 0.5*(vrtx2+vrtx3);
vrtx5 = 0.5*(vrtx3+vrtx1);
vrtx6 = 0.5*(vrtx1+vrtx2);
vrtx7 = (1/3)*(vrtx1+vrtx2+vrtx3);
% Short/small edges on control volume bndry. or inside primary triangles
NmlSmlEg = zeros(NumEms,3,2);  % Not unit-length; Thru rotating tan.vec.
NmlSmlEg(:,1,:) = (vrtx7-vrtx4) * [0,-1; 1,0];
NmlSmlEg(:,2,:) = (vrtx7-vrtx5) * [0,-1; 1,0];
NmlSmlEg(:,3,:) = (vrtx7-vrtx6) * [0,-1; 1,0];

%% 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);
% TO BE REVISED FOR EFFICIENCY 
DirichletNodeFlag = zeros(NumNds,1);
for ig=1:NumDirichletEgs
  k1 = TriMesh.edge(DirichletEdge(ig),1);
  k2 = TriMesh.edge(DirichletEdge(ig),2);
  DirichletNodeFlag(k1) = 1;
  DirichletNodeFlag(k2) = 1;
end
DirichletNode = find(DirichletNodeFlag);
NumDirichletNds = sum(DirichletNodeFlag);

%% Setting up 
DOFs = NumNds;
GlbMat = sparse(DOFs,DOFs);
GlbRHS = zeros(DOFs,1); 

%% Darcy_FVMD_TriP1: Assembling the global coefficient matrix 
% CP1BasGrad = zeros(NumEms,3,2);
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);
% For each elt., 3 bas.fxns., 3 small edges 
GradIntgrl = zeros(NumEms,3,3);
for i=1:3
  for j=1:3
    GradIntgrl(:,i,j) = dot(-[PermKs,PermKs].*squeeze(CP1BasGrad(:,i,:)),...
      squeeze(NmlSmlEg(:,j,:)),2);
  end
end
SignMat = [0 -1 1; 1 0 -1; -1 1 0];
EltStfMat = zeros(NumEms,3,3);
for i=1:3
  for j=1:3
    EltStfMat(:,i,j) = SignMat(i,1)*GradIntgrl(:,j,1)...
                     + SignMat(i,2)*GradIntgrl(:,j,2)...
                     + SignMat(i,3)*GradIntgrl(:,j,3);
  end
end
% Now assembly 
for i=1:3
  II = TriMesh.elem(:,i);
  for j=1:3
    JJ = TriMesh.elem(:,j);
    GlbMat = GlbMat + sparse(II,JJ,EltStfMat(:,i,j),DOFs,DOFs);
  end
end

%% Assembling the global right-hand side: GlbRHS 
% GlbRHS = zeros(DOFs,1);
EltRHS = zeros(NumEms,6);
% Applying the Gaussian quadrature on each small dual triangle 
NumQuadPts = size(GAUSSQUAD.TRIG,1);
% qp = zeros(NumEms,2);
% 
% For vertex #1: 1st triangle (1,6,7) 
x1 = vrtx1(:,1);  x2 = vrtx6(:,1);  x3 = vrtx7(:,1);
y1 = vrtx1(:,2);  y2 = vrtx6(:,2);  y3 = vrtx7(:,2);
area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * [x1 y1]...
     + GAUSSQUAD.TRIG(k,2) * [x2 y2]...
     + GAUSSQUAD.TRIG(k,3) * [x3 y3];
  EltRHS(:,1) = EltRHS(:,1) + GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp);
end
EltRHS(:,1) = EltRHS(:,1).*area;
% 
% For vertex #1: 2nd triangle (1,7,5) 
x1 = vrtx1(:,1);  x2 = vrtx7(:,1);  x3 = vrtx5(:,1);
y1 = vrtx1(:,2);  y2 = vrtx7(:,2);  y3 = vrtx5(:,2);
area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * [x1 y1]...
     + GAUSSQUAD.TRIG(k,2) * [x2 y2]...
     + GAUSSQUAD.TRIG(k,3) * [x3 y3];
  EltRHS(:,2) = EltRHS(:,2) + GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp);
end
EltRHS(:,2) = EltRHS(:,2).*area;
% 
% For vertex #2: 1st triangle (2,4,7) 
x1 = vrtx2(:,1);  x2 = vrtx4(:,1);  x3 = vrtx7(:,1);
y1 = vrtx2(:,2);  y2 = vrtx4(:,2);  y3 = vrtx7(:,2);
area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * [x1 y1]...
     + GAUSSQUAD.TRIG(k,2) * [x2 y2]...
     + GAUSSQUAD.TRIG(k,3) * [x3 y3];
  EltRHS(:,3) = EltRHS(:,3) + GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp);
end
EltRHS(:,3) = EltRHS(:,3).*area;
% 
% For vertex #2: 2nd triangle (2,7,6) 
x1 = vrtx2(:,1);  x2 = vrtx7(:,1);  x3 = vrtx6(:,1);
y1 = vrtx2(:,2);  y2 = vrtx7(:,2);  y3 = vrtx6(:,2);
area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * [x1 y1]...
     + GAUSSQUAD.TRIG(k,2) * [x2 y2]...
     + GAUSSQUAD.TRIG(k,3) * [x3 y3];
  EltRHS(:,4) = EltRHS(:,4) + GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp);
end
EltRHS(:,4) = EltRHS(:,4).*area;
%
% For vertex #3: 1st triangle (3,5,7) 
x1 = vrtx3(:,1);  x2 = vrtx5(:,1);  x3 = vrtx7(:,1);
y1 = vrtx3(:,2);  y2 = vrtx5(:,2);  y3 = vrtx7(:,2);
area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * [x1 y1]...
     + GAUSSQUAD.TRIG(k,2) * [x2 y2]...
     + GAUSSQUAD.TRIG(k,3) * [x3 y3];
  EltRHS(:,5) = EltRHS(:,5) + GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp);
end
EltRHS(:,5) = EltRHS(:,5).*area;
% 
% For vertex #3: 2nd triangle (3,7,4) 
x1 = vrtx3(:,1);  x2 = vrtx7(:,1);  x3 = vrtx4(:,1);
y1 = vrtx3(:,2);  y2 = vrtx7(:,2);  y3 = vrtx4(:,2);
area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * [x1 y1]...
     + GAUSSQUAD.TRIG(k,2) * [x2 y2]...
     + GAUSSQUAD.TRIG(k,3) * [x3 y3];
  EltRHS(:,6) = EltRHS(:,6) + GAUSSQUAD.TRIG(k,4)*EqnBC.fxnf(qp);
end
EltRHS(:,6) = EltRHS(:,6).*area;
% 
% Putting together 
for ie=1:NumEms
  for j=1:3
    GlbRHS(TriMesh.elem(ie,j)) = GlbRHS(TriMesh.elem(ie,j))...
      + EltRHS(ie,2*j-1) + EltRHS(ie,2*j);
  end
end

%% For reducing... 
% disp('Adjusting the global linear system...'); 
flag = zeros(DOFs,1);
flag(DirichletNode) = ones(NumDirichletNds,1);
FreeNd = find(~flag);
% sln vector to be used 
sln = zeros(DOFs,1);

%% Incorporating boundary conditions: Dirichlet as essential 
% Approach I: working on all these nodes simultaneously 
% Assuming NumDirichletNds>0
DiriNd = TriMesh.node(DirichletNode,:);
pDval = EqnBC.fxnpD(DiriNd);
sln(DirichletNode) = pDval;

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

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

return;