function sln = TransCDnv_IE_WG_TriP0P0RT0_AsmSlv(EqnBIC,CofRT0NB,...
  TriMesh,tyme,DirichletEdge,NeumannEdge,GAUSSQUAD)
%% Solving transient(transport) convection-diffusion eqn. (numer. velocity)
% by Implicit Euler + WG(P0;P0,RT0) on a triangular mesh
% James Liu, ColoState; 2012/07--2017/02 

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

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

%% Boundary info 
NumDirichletEgs = size(DirichletEdge,1);

%% Setup 
DOFs = NumEms + NumEgs;
% GlbMat = sparse(DOFs,DOFs);

%% Assembling matrices: For mass: 
GlbMatMass = sparse(DOFs,DOFs);
II = (1:NumEms);  JJ = II;
GlbMatMass = GlbMatMass + sparse(II,JJ,area,DOFs,DOFs);

%% Assembling matrices: For convection: 
GlbMatConv = sparse(DOFs,DOFs);
% Assuming numerical velocity in RT0 
GMRT0 = Hdiv_TriRT0_NmlzBas_GramMat(TriMesh);
CDWGB = WG_TriP0P0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(TriMesh);
AryConv = zeros(NumEms,4);
for i=1:4
  AryConv(:,i) = sum(CofRT0NB.*squeeze(CDWGB(:,i,:)).*GMRT0, 2);
end
% Assuming numerical velocity in BDM1: To be added later  
% Assuming numerical velocity in RT1: To be added later  
% Now assembly: Note the negative sign !!
II = (1:NumEms);  JJ = II;
GlbMatConv = GlbMatConv + sparse(II,JJ,-AryConv(:,1),DOFs,DOFs);
JJ = (1:NumEms);
for i=1:3
  II = NumEms + TriMesh.elem2edge(:,i);
  GlbMatConv = GlbMatConv + sparse(II,JJ,-AryConv(:,i+1),DOFs,DOFs);
end

%% Assembling matrices: For diffusion: 
GlbMatDiff = sparse(DOFs,DOFs);
D = EqnBIC.D;
DiffD(:,1,1) = D;
DiffD(:,2,2) = D;
GMK = Hdiv_TriRT0_NmlzBas_GramMatK(TriMesh,DiffD);
CDWGB = WG_TriP0P0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(TriMesh);
AryEE = zeros(NumEms,1);
AryEG = zeros(NumEms,3);
AryGG = zeros(NumEms,3,3);
% JL20151225: TO BE REVISED FOR EFFECIENCY !!
for ie=1:NumEms
  EltGM = squeeze(GMK(ie,:,:));
  AryEE(ie) = squeeze(CDWGB(ie,1,:))' * EltGM * squeeze(CDWGB(ie,1,:));
  for i=1:3
    AryEG(ie,i) = squeeze(CDWGB(ie,1,:))' * EltGM * squeeze(CDWGB(ie,1+i,:));
    for j=1:3
      AryGG(ie,i,j) = squeeze(CDWGB(ie,1+i,:))' * EltGM * squeeze(CDWGB(ie,1+j,:));
    end
  end
end
% Part 1: element-element interaction 
II = (1:NumEms);  JJ = II;
GlbMatDiff = GlbMatDiff + sparse(II,JJ,AryEE,DOFs,DOFs);
% Part 2: element-edge interaction 
II = (1:NumEms);
for j=1:3
  JJ = NumEms + TriMesh.elem2edge(:,j);
  GlbMatDiff = GlbMatDiff + sparse(II,JJ,AryEG(:,j),DOFs,DOFs);
  GlbMatDiff = GlbMatDiff + sparse(JJ,II,AryEG(:,j),DOFs,DOFs);
end
% Part 3: edge-edge interaction 
for i=1:3
  II = NumEms + TriMesh.elem2edge(:,i);
  for j=i:3  % Utilizing symmetry 
    JJ = NumEms + TriMesh.elem2edge(:,j);
    if (j==i) 
      GlbMatDiff = GlbMatDiff + sparse(II,JJ,AryGG(:,i,i),DOFs,DOFs);
    else
      GlbMatDiff = GlbMatDiff + sparse(II,JJ,AryGG(:,i,j),DOFs,DOFs);
      GlbMatDiff = GlbMatDiff + sparse(JJ,II,AryGG(:,i,j),DOFs,DOFs);
    end
  end
end
 
%% Assembling matrices: Combined 
GlbMat = GlbMatMass + Deltat*(GlbMatConv+GlbMatDiff);
clear  GlbMatConv  GlbMatDiff;

%% Initial numerical solution 
sold = zeros(DOFs,1);
% For element interior 
sold(1:NumEms) = zeros(NumEms,1);
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),:);
  sold(1:NumEms) = sold(1:NumEms) ...
    + GAUSSQUAD.TRIG(k,4) * EqnBIC.fxnc0(qp);  % */area 
end
% For mesh skeleton 
x1 = TriMesh.node(TriMesh.edge(:,1),1);
y1 = TriMesh.node(TriMesh.edge(:,1),2);
x2 = TriMesh.node(TriMesh.edge(:,2),1);
y2 = TriMesh.node(TriMesh.edge(:,2),2);
NumQuadPts = size(GAUSSQUAD.LINE,1);
qp = zeros(NumEgs,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;
  sold(NumEms+(1:NumEgs)) = sold(NumEms+(1:NumEgs)) ... 
    + GAUSSQUAD.LINE(k,3)*EqnBIC.fxnc0(qp);  % */ LenEg
end 
% 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,1);
  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),:);
    EltSrc = EltSrc + GAUSSQUAD.TRIG(k,4) * EqnBIC.fxnf(qp,tnew);
  end
  EltSrc = Deltat * (area.*EltSrc) ;
  GlbRHS(1:NumEms,1) = GlbRHS(1:NumEms,1) + EltSrc;

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

  % For reducing ... 
  flag = zeros(DOFs,1);
  flag(NumEms+DirichletEdge) = ones(NumDirichletEgs,1);
  EmFreeEg = find(~flag);
  snew = zeros(DOFs,1);

  % Incorporating boundary conditions: Dirichlet as essential 
  % Approach I: working on all these edges simultaneously 
  if NumDirichletEgs>0 
    cDavg = zeros(NumDirichletEgs,1);
    x1 = TriMesh.node(TriMesh.edge(DirichletEdge,1),1);
    y1 = TriMesh.node(TriMesh.edge(DirichletEdge,1),2);
    x2 = TriMesh.node(TriMesh.edge(DirichletEdge,2),1);
    y2 = TriMesh.node(TriMesh.edge(DirichletEdge,2),2);
    NumQuadPts = size(GAUSSQUAD.LINE,1);
    qp = zeros(NumDirichletEgs,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;
      cDavg = cDavg + GAUSSQUAD.LINE(k,3)*EqnBIC.fxncD(qp,tnew);
    end
    snew(NumEms+DirichletEdge) = cDavg;
  end

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

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

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

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

return;