function  sln = Parab_IE_CG_TriP1_AsmSlv(...
  EqnBIC,TriMesh,tyme,DirichletEdge,NeumannEdge,GAUSSQUAD)
%% Solving the parabolic equation 
% by Implicit Euler + Continuous Galerkin on a triangular mesh
% James Liu, ColoState; 2012/07--2017/02 

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

%% Tyme info
NT = tyme.NT;  MT = tyme.MT;  Deltat = tyme.Deltat;

%% 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,NT+1);

%% Assembling the "static" global matrix: Mass 
% For mass: elementwise 
% Alternative: EltMatMass = CG_TriP1_GramMat(TriMesh);
EltMatMass = zeros(NumEms,3,3);
for i=1:3
  for j=1:3
    if (j==i) 
      EltMatMass(:,i,j) = area/6;
    else 
      EltMatMass(:,i,j) = area/12;
    end
  end
end
% For mass: global 
GlbMatMass = sparse(DOFs,DOFs);
for i=1:3
  II = TriMesh.elem(:,i);
  for j=1:3
    JJ = TriMesh.elem(:,j);
    GlbMatMass = GlbMatMass + sparse(II,JJ,EltMatMass(:,i,j),DOFs,DOFs);
  end
end

%% Assembling the "static" global matrix: Stiffness/Diffusion  
% For diffusion: elementwise 
% Alternative: EltMatDiff = CG_TriP1_GradGradMat(TriMesh); 
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);
P1BasGrad = zeros(NumEms,3,2);
P1BasGrad(:,1,:) = [(y2-y3)./(2*area), (x3-x2)./(2*area)];
P1BasGrad(:,2,:) = [(y3-y1)./(2*area), (x1-x3)./(2*area)];
P1BasGrad(:,3,:) = [(y1-y2)./(2*area), (x2-x1)./(2*area)];
EltMatDiff = zeros(NumEms,3,3);
for i=1:3
  for j=1:3
    EltMatDiff(:,i,j) = squeeze(dot(P1BasGrad(:,i,:),P1BasGrad(:,j,:),3))...
      .* area;
  end
end
% For diffusion: global 
GlbMatDiff = sparse(DOFs,DOFs);
for i=1:3
  II = TriMesh.elem(:,i);
  for j=i:3  % Utilizing symmetry
    JJ = TriMesh.elem(:,j);
    if (j==i) 
      GlbMatDiff = GlbMatDiff + sparse(II,JJ,EltMatDiff(:,i,j),DOFs,DOFs);
    else
      GlbMatDiff = GlbMatDiff + sparse(II,JJ,EltMatDiff(:,i,j),DOFs,DOFs);
      GlbMatDiff = GlbMatDiff + sparse(JJ,II,EltMatDiff(:,i,j),DOFs,DOFs);
    end
  end
end

%% Assembling the "static" global matrix: Combined 
GlbMat = GlbMatMass + Deltat * GlbMatDiff;
% clear 

%% Initial numerical solution
sold = zeros(DOFs,1);
% Setting the initial numerical solution by interpolation 
sold = EqnBIC.fxnc0(TriMesh.node);
% Storing 
sln(:,1) = sold;

%% Time-marching 
for n=1:MT
  tnew = tyme.tm(n+1);

  % For the global right-hand side 
  GlbRHS = zeros(DOFs,1);

  % Assembling the source term at the new time 
  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),:);
    % Note: Lagrange P1 bas.fxn.val. = barycentric.crd. = Quad.weight 
    for j=1:3
      EltSrc(:,j) = EltSrc(:,j) + GAUSSQUAD.TRIG(k,4)...
        * GAUSSQUAD.TRIG(k,j).*EqnBIC.fxnf(qp,tnew);
    end
  end
  EltSrc = EltSrc .* ([area,area,area]*Deltat);
  % JL20160321: TO BE REVISED FOR EFFECIENCY 
  % Putting together 
  for ie=1:NumEms
    for j=1:3
      GlbRHS(TriMesh.elem(ie,j)) = ...
        GlbRHS(TriMesh.elem(ie,j)) + EltSrc(ie,j);
    end
  end

  % Assembling the old mass 
  GlbRHS = GlbRHS + GlbMatMass*sold;

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

  % Incorporating boundary conditions: Dirichlet as essential 
  % Approach I: working on all these nodes simultaneously 
  % Assuming NumDirichletNds>0
  DiriNd = TriMesh.node(DirichletNode,:);
  cDval = EqnBIC.fxncD(DiriNd,tnew);
  snew(DirichletNode) = cDval;

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

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

  % Storing  
  sln(:,n+1) = snew;

  % Marching to the next time step  
  sold = snew;
end

return;