function sln = SteadyCD_CG_TriP1_AsmSlv(...
  EqnBC,TriMesh,CofRT0EgBas,NmlFlux,DirichletEdge,NeumannEdge,GAUSSQUAD) 
%% Solving steady convection-diffusion eqn. by CG P1 on a triangular mesh 
% James Liu, ColoState; 2012/07--2017/02

%% EqnBC
D = EqnBC.D; 

%% Mesh 
NumNds = TriMesh.NumNds;
NumEms = TriMesh.NumEms;
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);
  DirichletNodeFlag(k1) = 1;
  k2 = TriMesh.edge(DirichletEdge(ig),2);
  DirichletNodeFlag(k2) = 1;
end
DirichletNode = find(DirichletNodeFlag);
NumDirichletNds = sum(DirichletNodeFlag);

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

%% Assembling the global matrices 
% For diffusion 
% CP1BasGrad = zeros(NumEms,3,2);
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);
EltMat = zeros(NumEms,3,3);
for i=1:3
  for j=1:3
    EltMat(:,i,j) = dot(squeeze(CP1BasGrad(:,i,:)),...
      squeeze(CP1BasGrad(:,j,:)),2) .* area * D;
  end
end
GlbMatDiff = sparse(NumNds,NumNds);
for i=1:3
  II = TriMesh.elem(:,i);
  for j=1:3
    JJ = TriMesh.elem(:,j);
    GlbMatDiff = GlbMatDiff + sparse(II,JJ,EltMat(:,i,j),NumNds,NumNds);
  end
end

% For the convection 
GlbMatConv = CG_TriP1_GlbMatConv(TriMesh,CofRT0EgBas);
% % For the flux 
% GMF = CG_TriP1_GlbMatFlux(TriMesh,NmlFlux);
% All put together 
GlbMat = GlbMatDiff - GlbMatConv;

%% Assembling the global right-hand side 
% Applying the Gaussian quadrature on triangular elements 
EltSrc = zeros(NumEms,3);
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),:);
  fval = EqnBC.fxnf(qp);
  % Note: The Lagrange P1 bas.fxn.vals. = Gauss quadrature weights
  for j=1:3
    EltSrc(:,j) = EltSrc(:,j)...
      + GAUSSQUAD.TRIG(k,4) * fval .* GAUSSQUAD.TRIG(k,1);
  end
end
EltSrc = EltSrc.*[area,area,area];
% JL20160326: NEED TO FIND A MORE EFFECIENT WAY!
for ie=1:NumEms
  for j=1:3
    GlbRHS(TriMesh.elem(ie,j)) = GlbRHS(TriMesh.elem(ie,j)) + EltSrc(ie,j);
  end
end

%% For reducing... 
flag = zeros(DOFs,1);
flag(DirichletNode) = ones(NumDirichletNds,1);
FreeNd = find(~flag);

%% Incorporating boundary conditions: Dirichlet as essential 
% Assuming NumDirichletNds>0
DiriNd = TriMesh.node(DirichletNode,:);
cDval = EqnBC.fxncD(DiriNd);
sln(DirichletNode) = cDval;

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

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

return;