function sln = TransCDRnv_ELLAM_CG_RectQ1_AsmSlv(...
  EqnBIC,CofRT0NB,RectMesh,tyme,DirichletEdge,NeumannEdge,GAUSSQUAD)
%% Solving transient convection-diffusion-reaction eqn. (numer. velocity)
% by ELLAM + CG.Q1 on a rectangular mesh
% Assuming a constant reaction coefficient
% James Liu, ColoState;  2012/07--2017/02 

%% EqnBIC
D = EqnBIC.D;
R = EqnBIC.R;

%% Mesh info 
NumNds = RectMesh.NumNds;
NumEms = RectMesh.NumEms;
% NumEgs = TriMesh.NumEgs;
% area = TriMesh.area;
hx = RectMesh.hx;  hy = RectMesh.hy;
xa = RectMesh.xa;  % xb = RectMesh.xb;  
yc = RectMesh.yc;  % yd = RectMesh.yd;

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

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

%% Setup 
DOFs = NumNds;
GlbMat = sparse(DOFs,DOFs);
GlbRHS = zeros(DOFs,1);
sln = zeros(DOFs,NT+1);

%% Assembling "static" global matrices 
EGM = CG_RectQ1_GramMat(RectMesh);  % For mass 
EGGM = CG_RectQ1_GradGradMat(RectMesh);  % For diffusion 
for i=1:4
  II = RectMesh.elem(:,i);
  for j=1:4
    JJ = RectMesh.elem(:,j);
    Mij = EGM(:,i,j);
    GlbMat = GlbMat + sparse(II,JJ,Mij,NumNds,NumNds);
    Sij = Deltat * D * EGGM(:,i,j);  % Aprx.: Deltat indep. of quad.pts. 
    GlbMat = GlbMat + sparse(II,JJ,Sij,NumNds,NumNds);
  end
end
clear  EGM  EGGM  Mij  Sij;

%% Initial numerical solution 
sold = CG_RectQ1_intpln(EqnBIC.fxnc0, RectMesh);

%% Storing 
sln(:,1) = sold;

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

  % GlbRHS: There are three big parts 
  GlbRHS = zeros(NumNds,1);

  % GlbRHS Part I: Computing and assembling the old mass 
  EltOldMass = zeros(NumEms,4);  % 4 Lagrangian Q1 bas.fxns.
  % 
  NumQuadPts = size(GAUSSQUAD.RECT,1);
  GAUSSQUADRECT1 = ones(size(GAUSSQUAD.RECT)) - GAUSSQUAD.RECT;
  qp = zeros(NumEms,2);
  for k=1:NumQuadPts 
    qp(:,1) = GAUSSQUAD.RECT(k,1) * RectMesh.node(RectMesh.elem(:,1),1)...
            + GAUSSQUADRECT1(k,1) * RectMesh.node(RectMesh.elem(:,3),1);
    qp(:,2) = GAUSSQUAD.RECT(k,2) * RectMesh.node(RectMesh.elem(:,1),2)...
            + GAUSSQUADRECT1(k,2) * RectMesh.node(RectMesh.elem(:,3),2);
    % Non-normalized 1-dim Lagrangian bas.fxns., See also comments later 
    X1 = RectMesh.node(RectMesh.elem(:,3),1) - qp(:,1);
    Y1 = RectMesh.node(RectMesh.elem(:,3),2) - qp(:,2);
    X2 = qp(:,1) - RectMesh.node(RectMesh.elem(:,1),1);
    Y2 = qp(:,2) - RectMesh.node(RectMesh.elem(:,1),2);
    % 
    % JL20151212: TO BE INVESTIGATED FURTHER 
    % Back-tracking the quadrature points: Euler 
    X = qp(:,1) - RectMesh.EmCntr(:,1);
    Y = qp(:,2) - RectMesh.EmCntr(:,2);
    vel(:,1) = CofRT0NB(:,1) + CofRT0NB(:,3).*X;
    vel(:,2) = CofRT0NB(:,2) + CofRT0NB(:,4).*Y;
    foot = qp - Deltat * vel;
    % 
    % JL20151212: TO BE REVISED 
    % Locating the foot under assumption of uniform mesh, w/ modification 
    Istar = floor((foot(:,1)-xa)/hx) + 1;
    Jstar = floor((foot(:,2)-yc)/hy) + 1;
    kk = find(Istar<1);  Istar(kk) = 1;
    kk = find(Istar>RectMesh.nx);  Istar(kk) = RectMesh.nx;
    kk = find(Jstar<1);  Jstar(kk) = 1;
    kk = find(Jstar>RectMesh.ny);  Jstar(kk) = RectMesh.ny;
    [Istar, Jstar];
    Kstar = (Istar-1)*RectMesh.ny + Jstar;
    % 
    % Lagrangian Q1-interpolatiion 
    k1 = RectMesh.elem(Kstar,1);  k2 = RectMesh.elem(Kstar,2);  
    k3 = RectMesh.elem(Kstar,3);  k4 = RectMesh.elem(Kstar,4);
    x1 = RectMesh.node(k1,1);  y1 = RectMesh.node(k1,2);
    x2 = RectMesh.node(k3,1);  y2 = RectMesh.node(k3,2);
    Xstar = (foot(:,1)-x1)./(x2-x1);  
    Ystar = (foot(:,2)-y1)./(y2-y1);
    Cstar = sold(k1) .* (1-Xstar).*(1-Ystar) ...
          + sold(k2) .*     Xstar.*(1-Ystar) ... 
          + sold(k3) .*     Xstar.*Ystar ... 
          + sold(k4) .* (1-Xstar).*Ystar ;
    % 
    % Incorporating the constant reaction coefficient
    Cstar = Cstar * exp(-R*Deltat);
    % 
    % JL20151212: TO BE ADDED LATER 
    % Computing Jacobian J(foot,qp) 
    % 
    % Integration: Element area is not needed due to the form of bas.fxns.
    EltOldMass(:,1) = EltOldMass(:,1) + GAUSSQUAD.RECT(k,3) * (X1.*Y1).*Cstar;
    EltOldMass(:,2) = EltOldMass(:,2) + GAUSSQUAD.RECT(k,3) * (X2.*Y1).*Cstar;
    EltOldMass(:,3) = EltOldMass(:,3) + GAUSSQUAD.RECT(k,3) * (X2.*Y2).*Cstar;
    EltOldMass(:,4) = EltOldMass(:,4) + GAUSSQUAD.RECT(k,3) * (X1.*Y2).*Cstar;
  end
  % Assembling 
  for j=1:4
    i = RectMesh.elem(:,j);
    tmp = GlbRHS(i);
    GlbRHS(i) = tmp + EltOldMass(:,j);
  end
  
  % GlbRHS Part II: Computing and assembling elementwise source 
  EltSrcNew = zeros(NumEms,4);
  NumQuadPts = size(GAUSSQUAD.RECT,1);
  GAUSSQUADRECT1 = ones(size(GAUSSQUAD.RECT)) - GAUSSQUAD.RECT;
  qp = zeros(NumEms,2);
  for k=1:NumQuadPts 
    qp(:,1) = GAUSSQUAD.RECT(k,1) * RectMesh.node(RectMesh.elem(:,1),1)...
            + GAUSSQUADRECT1(k,1) * RectMesh.node(RectMesh.elem(:,3),1);
    qp(:,2) = GAUSSQUAD.RECT(k,2) * RectMesh.node(RectMesh.elem(:,1),2)...
            + GAUSSQUADRECT1(k,2) * RectMesh.node(RectMesh.elem(:,3),2);
    fval = EqnBIC.fxnf(qp,tnew);
    % Non-normalized 1-dim Lagrangian bas.fxns., See also comments later 
    X1 = RectMesh.node(RectMesh.elem(:,3),1) - qp(:,1);
    Y1 = RectMesh.node(RectMesh.elem(:,3),2) - qp(:,2);
    X2 = qp(:,1) - RectMesh.node(RectMesh.elem(:,1),1);
    Y2 = qp(:,2) - RectMesh.node(RectMesh.elem(:,1),2);
    EltSrcNew(:,1) = EltSrcNew(:,1) + GAUSSQUAD.RECT(k,3)*fval.*(X1.*Y1);
    EltSrcNew(:,2) = EltSrcNew(:,2) + GAUSSQUAD.RECT(k,3)*fval.*(X2.*Y1);
    EltSrcNew(:,3) = EltSrcNew(:,3) + GAUSSQUAD.RECT(k,3)*fval.*(X2.*Y2);
    EltSrcNew(:,4) = EltSrcNew(:,4) + GAUSSQUAD.RECT(k,3)*fval.*(X1.*Y2);
  end
  % Assembling 
  for j=1:4
    i = RectMesh.elem(:,j);
    tmp = GlbRHS(i);
    GlbRHS(i) = tmp + Deltat * EltSrcNew(:,j);
  end

  % GlbRHS Part III: Incorporating boundary conditions: Neumann 
  % None for now 

  % Finding Dirichlet nodes 
  % JL20151212: TO BE REVISED FOR EFFICIENCY !
  flag = zeros(NumNds,1);
  for ig=1:NumDirichletEgs
    k1 = RectMesh.edge(DirichletEdge(ig),1);
    k2 = RectMesh.edge(DirichletEdge(ig),2);
    flag(k1) = 1;
    flag(k2) = 1;
  end
  DirichletNode = find(flag);
  FreeNd = find(~flag);
  snew = zeros(NumNds,1);

  % Incorporating boundary conditions: Dirichlet as essential 
  cDval = EqnBIC.fxncD(RectMesh.node(DirichletNode,:),tnew);
  snew(DirichletNode) = cDval;

  % Reducing 
  GlbRHS = GlbRHS - GlbMat*snew;

  %% 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;