function [sln] = Darcy_CG_TriP1_AsmSlv(...
  EqnBC,TriMesh,PermK,DirichletEdge,NeumannEdge,GAUSSQUAD)
%% Darcy: CG: P1 on a triangular mesh: Assembling & Solving 
% James Liu, ColoState; 2012/07--2018/12 

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

%% 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 the global linear system 
DOFs = NumNds;
% GlbMat = sparse(DOFs,DOFs);
% GlbRHS = zeros(DOFs,1);
% sln = zeros(DOFs,1);

%% Assembling the global coefficient matrix: GlbMat 
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);
GlbMat = sparse(DOFs,DOFs);
for i=1:3
  II = TriMesh.elem(:,i);
  for j=1:3
    JJ = TriMesh.elem(:,j);
    Aij = CP1BasGrad(:,i,1) .* PermK(:,1,1) .* CP1BasGrad(:,j,1) ... 
        + CP1BasGrad(:,i,1) .* PermK(:,1,2) .* CP1BasGrad(:,j,2) ...
        + CP1BasGrad(:,i,2) .* PermK(:,2,1) .* CP1BasGrad(:,j,1) ...
        + CP1BasGrad(:,i,2) .* PermK(:,2,2) .* CP1BasGrad(:,j,2) ;
    Aij = Aij.*area;
    GlbMat = GlbMat + sparse(II,JJ,Aij,NumNds,NumNds);
  end
end

% %% JL20131004: NEED DOUBLE-CHECK
% %% 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),:);
%   BaryCrd = TriMesh_BaryCrd(TriMesh,qp);
%   fval = EqnBC.fxnf(qp);
%   for j=1:3
%     GlbRHS(TriMesh.elem(:,j)) = GlbRHS(TriMesh.elem(:,j))...
%       + GAUSSQUAD.TRIG(k,4) * (fval.*BaryCrd(:,j)) .* area;
%   end  
% end

%% JL20131004: Old style temporarily 
%% Assembling GlbRHS 
GlbRHS = zeros(DOFs,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
phi = zeros(1,3);
% A = zeros(1,3);
for ie=1:NumEms
  % k1 = TriMesh.elem(ie,1);  k2 = TriMesh.elem(ie,2);  k3 = TriMesh.elem(ie,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);
  % AA = 0.5*((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
  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),:);
    fval = EqnBC.fxnf(qp);
    % xx = qp(:,1);  yy = qp(:,2);
    % A(1) = 0.5*((x2-xx).*(y3-yy)-(x3-xx).*(y2-yy));
    % A(2) = 0.5*((xx-x1).*(y3-y1)-(x3-x1).*(yy-y1));
    % A(3) = 0.5*((x2-x1).*(yy-y1)-(xx-x1).*(y2-y1));
    % for i=1:3;  phi(i) = A(i)/AA;  end  
    for j=1:3;  phi(j) = GAUSSQUAD.TRIG(k,j);  end
    for j=1:3
      GlbRHS(TriMesh.elem(ie,j)) = GlbRHS(TriMesh.elem(ie,j))...
        + GAUSSQUAD.TRIG(k,4) * (fval*phi(j)) * area(ie);
    end
  end
end

%% Incorporating boundary conditions: Neumann as natural 
% Approach I: working on all these edges simultaneously 
if NumNeumannEgs>0
  k1 = TriMesh.edge(NeumannEdge,1);  
  k2 = TriMesh.edge(NeumannEdge,2);
  x1 = TriMesh.node(k1,1);  y1 = TriMesh.node(k1,2);
  x2 = TriMesh.node(k2,1);  y2 = TriMesh.node(k2,2);
  % LenEgNeum = 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;
    phi1 = sqrt((x2-qp(:,1)).^2+(y2-qp(:,2)).^2);  % ./LenEgNeum;
    phi2 = sqrt((qp(:,1)-x1).^2+(qp(:,2)-y1).^2);  % ./LenEgNeum;
    uNval = EqnBC.fxnuN(qp);
    GlbRHS(k1) = GlbRHS(k1) - GAUSSQUAD.LINE(k,3)*(uNval.*phi1);
    GlbRHS(k2) = GlbRHS(k2) - GAUSSQUAD.LINE(k,3)*(uNval.*phi2);
  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;