function [NumerDisp,NumerPres,NumerVelCof] ...
  = LinPoroElas_MS_TriIE_CGEP1_WGP0P0RT0_AsmSlv( ...
    EqnBIC,BndryCondTypeElas,BndryCondTypeDarcy,TriMesh,tyme,GAUSSQUAD)
%% Solving linear poroelasticity as a monolithic system (MS) 
% on a triangular mesh by implicit Euler with 
% -- CG Bernardi-Raugel(1) for linear elasticity (solid displacement)
% -- WG(P0,P0;RT0) for Darcy flow (fluid pressure) 
% Zhuoran Wang, SYSU; James Liu ColoState; 2020/09--2021/04 

%% Equation info 
lambda = EqnBIC.lambda;
    mu = EqnBIC.mu;
 alpha = EqnBIC.alpha;
    c0 = EqnBIC.c0;

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

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

%% 
% For Elas: Sorting out boundary edges: Dirichlet, Neumann 
DirichletEdgeElas = find(BndryCondTypeElas(TriMesh.BndryEdge+1)==1);
NeumannEdgeElas   = find(BndryCondTypeElas(TriMesh.BndryEdge+1)==2);

NumDirichletEgsElas = size(DirichletEdgeElas,1);
NumNeumannEgsElas = size(NeumannEdgeElas,1);


% JL20160808: TO BE REVISED FOR EFFICIENCY !
DirichletNodeElasFlag = zeros(NumNds,1);
for ig=1:NumDirichletEgsElas 
  k1 = TriMesh.edge(DirichletEdgeElas(ig),1);
  k2 = TriMesh.edge(DirichletEdgeElas(ig),2);
  DirichletNodeElasFlag(k1) = 1;
  DirichletNodeElasFlag(k2) = 1;
end
DirichletNodeElas = find(DirichletNodeElasFlag);
NumDirichletNdsElas = sum(DirichletNodeElasFlag);

% For Darcy: Sorting out boundary edges: Dirichlet, Neumann 
DirichletEdgeDarcy = find(BndryCondTypeDarcy(TriMesh.BndryEdge+1)==1);
NeumannEdgeDarcy   = find(BndryCondTypeDarcy(TriMesh.BndryEdge+1)==2);
NumDirichletEgsDarcy = size(DirichletEdgeDarcy,1);
NumNeumannEgsDarcy = size(NeumannEdgeDarcy,1);

%% Setting up the discrete linear system 
DOFsDspl = 2*NumNds + NumEgs;
DOFsPres = NumEms + NumEgs;
NumerDisp = zeros(DOFsDspl,NT+1);
NumerPres = zeros(DOFsPres,NT+1);
NumerVelCof = zeros(NumEms,4,NT);

%% Assembling GlbMatElas 
% Element matrix: strain-strain 
EltMatStrnStrn = zeros(NumEms,9,9);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
%   xhat = GAUSSQUADRECT1(k,1);
%   yhat = GAUSSQUADRECT1(k,2);
%   jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
%       - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh,...
    GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  StrnPhi = zeros(NumEms,9,2,2);
  StrnPhi(:,:,1,1) = BFGrad(:,:,1,1);
  StrnPhi(:,:,2,2) = BFGrad(:,:,2,2);
  StrnPhi(:,:,1,2) = 0.5*(BFGrad(:,:,1,2)+BFGrad(:,:,2,1));
  StrnPhi(:,:,2,1) = StrnPhi(:,:,1,2);
  for i=1:9
    for j=1:9
      EltMatStrnStrn(:,i,j) = EltMatStrnStrn(:,i,j)...
        + ( StrnPhi(:,i,1,1).*StrnPhi(:,j,1,1)...
          + StrnPhi(:,i,1,2).*StrnPhi(:,j,1,2)...
          + StrnPhi(:,i,2,1).*StrnPhi(:,j,2,1)...
          + StrnPhi(:,i,2,2).*StrnPhi(:,j,2,2)) .* area * GAUSSQUAD.TRIG(k,4);
    end
  end
end
% Element matrix: div(avg)-div(avg) 
EltMatDivDiv = zeros(NumEms,9,9);
EltDivAvg = zeros(NumEms,9);
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);
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),:);
  [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh, ...
      GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  for j=1:9
    EltDivAvg(:,j) = EltDivAvg(:,j) ...
      +  BFDiv(:,j) .*  GAUSSQUAD.TRIG(k,4);
  end
end
for i=1:9
  for j=1:9
    EltMatDivDiv(:,i,j) = EltDivAvg(:,i) .* EltDivAvg(:,j) .* area;
  end
end
% Altogether 
% EltMat = zeros(NumEms,9,9);
EltMat = (2*mu)*EltMatStrnStrn + lambda*EltMatDivDiv;
% Modification for signs of edge normals 
for i=1:3
  for j=1:9
    EltMat(:,6+i,j) = EltMat(:,6+i,j) .* TriMesh.EmEgSign(:,i);
  end
end
for i=1:9
  for j=1:3
    EltMat(:,i,6+j) = EltMat(:,i,6+j) .* TriMesh.EmEgSign(:,j);
  end
end
% Now assembling GlbMatElas 
GlbMatElas = sparse(DOFsDspl,DOFsDspl);
pos = zeros(NumEms,9);
for j=1:3  % Treating vertices and edges in one loop 
  k = TriMesh.elem(:,j);
  pos(:,2*j-1) = 2*k - 1;
  pos(:,2*j  ) = 2*k;
  k = TriMesh.elem2edge(:,j);
  pos(:,6+j) = 2*NumNds + k;
end
for i=1:9 
  II = pos(:,i);
  for j=1:9 
    JJ = pos(:,j);
    GlbMatElas = GlbMatElas + sparse(II,JJ, EltMat(:,i,j), DOFsDspl,DOFsDspl);
  end
end

%% Assembling GlbMatDarcy 
% GlbMatDarcyMass = sparse(DOFsPres,DOFsPres);
II = (1:NumEms);  JJ = II;
GlbMatDarcyMass = sparse(II,JJ, area, DOFsPres,DOFsPres);
% Sampling permeability 
PermK = Darcy_SmplnPerm_TriMesh(EqnBIC.fxnK, TriMesh, GAUSSQUAD);
GMK = Hdiv_TriRT0_NmlzBas_GramMatK(TriMesh,PermK);
CDWGB = WG_TriP0P0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(TriMesh);
EltMat = zeros(NumEms,4,4);
% JL20210414: TO BE REVISED FOR EFFECIENCY 
for ie=1:NumEms
  EltGM = squeeze(GMK(ie,:,:));
  for i=1:4
    for j=1:4
      EltMat(ie,i,j) = squeeze(CDWGB(ie,i,:))' *EltGM* squeeze(CDWGB(ie,j,:));
    end
  end
end
pos = zeros(NumEms,4);
pos(:,1) = (1:NumEms);
for j=1:3
  pos(:,1+j) = NumEms + TriMesh.elem2edge(:,j);
end
GlbMatDarcyStf = sparse(DOFsPres,DOFsPres);
for i=1:4
  II = pos(:,i);
  for j=1:4
    JJ = pos(:,j);
    GlbMatDarcyStf = GlbMatDarcyStf ...
      + sparse(II,JJ, EltMat(:,i,j), DOFsPres,DOFsPres);
  end
end
% Altogether 
GlbMatDarcy = c0*GlbMatDarcyMass + Deltat*GlbMatDarcyStf;
clear EltMat;

%% JL20171231: TO BE REVISED 
%% Assembling GlbMatDarcyElas 
% Coupling: P0 (test fxns.) vs. averages of BR1 (trial fxns.) 
EltMatP0DivBR1 = zeros(NumEms,9);
for j=1:9
  EltMatP0DivBR1(:,j) = EltDivAvg(:,j) .* area;
end
% Adjusting for signs of edge normal vectors 
for j=1:3
  EltMatP0DivBR1(:,6+j) = EltMatP0DivBR1(:,6+j) .* TriMesh.EmEgSign(:,j);
end
GlbMatDarcyElas = sparse(DOFsPres,DOFsDspl);
pos = zeros(NumEms,9);
for j=1:3  % For nodes and edges (treated simultaneously) 
 pos(:,2*j-1) = 2*TriMesh.elem(:,j) - 1;
 pos(:,2*j)   = 2*TriMesh.elem(:,j);
 pos(:,6+j) = 2*NumNds + TriMesh.elem2edge(:,j);
end
II = (1:NumEms);
for j=1:9
  JJ = pos(:,j);
  GlbMatDarcyElas = GlbMatDarcyElas ... 
    + sparse(II, JJ, alpha*EltMatP0DivBR1(:,j), DOFsPres, DOFsDspl);
end

%% Producing GlbMatElasDarcy: Coupling Elas (test fxns.) w/ Darcy (trial fxns.)
GlbMatElasDarcy = -GlbMatDarcyElas';

%% MS: The monolithic system 
GlbMat = [GlbMatElas,  GlbMatElasDarcy; ...
          GlbMatDarcyElas, GlbMatDarcy];

%% Computing projections/interpolations 
% in mandel, just 0 as initial conditions.
% [ProjPres,IntpltDspl] = LinPoroElas_TriMesh_ProjPresP0P0_IntpltDsplBR1(...
%   EqnBIC,TriMesh,tyme,GAUSSQUAD);
% Initial conditions: 
% NumerPres(:,1) = ProjPres(:,1);
% NumerDspl(:,1) = IntpltDspl(:,1);
NumerPres(:,1) = 0;
NumerDisp(:,1) = 0;

%% Now time-marching 
fprintf('Time-marching step: ');
StepInTimeMarching = 0;
for n=1:MT
  fprintf(repmat('\b', 1, StepInTimeMarching));
  StepInTimeMarching = fprintf('%d', n);
  tnew = tyme.tm(n+1);

  % Assembling GlbRHS: For solid 
  % Computing elementwise contributions of 12 basis functions 
  EltRHS = zeros(NumEms,9);
  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 = EqnBIC.fxnf(qp,tnew);
    BFVal = CG_TriBR1_BasFxnVal(TriMesh,...
      GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
    for j=1:9
      EltRHS(:,j) = EltRHS(:,j) + GAUSSQUAD.TRIG(k,4)...
        *( BFVal(:,j,1).*fval(:,1) + BFVal(:,j,2).*fval(:,2) ).*area;
    end
  end
  % Modification for 3 edge-based basis functions 
  for j=1:3
    EltRHS(:,6+j) = EltRHS(:,6+j) .* TriMesh.EmEgSign(:,j);
  end

  GlbRhsDspl = zeros(DOFsDspl,1);
  for ie=1:NumEms
    for j=1:3  % Treating vertices and edges in one loop 
      k = TriMesh.elem(ie,j);
      GlbRhsDspl(2*k-1) = GlbRhsDspl(2*k-1) + EltRHS(ie,2*j-1);
      GlbRhsDspl(2*k  ) = GlbRhsDspl(2*k  ) + EltRHS(ie,2*j  );
      k = TriMesh.elem2edge(ie,j);
      GlbRhsDspl(2*NumNds+k) = GlbRhsDspl(2*NumNds+k) + EltRHS(ie,6+j);
    end
  end
  
  % Assembling GlbRHS: For fluid 
  GlbRhsPres = zeros(DOFsPres,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),:);
    sval = EqnBIC.fxns(qp,tnew);
    GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms)...
      + GAUSSQUAD.TRIG(k,4) * sval.* area * Deltat;
  end
  
  % Assembling GlbRHS: Pressure at old time 
  GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) ...
    + c0 * area.*squeeze(NumerPres(1:NumEms,n));

  % Assembling GlbRHS: Solid displacement divergence term at old time 
  tmp = zeros(NumEms,9);
  for j=1:3  % Treating vertices and edges in one loop 
    k = TriMesh.elem(:,j);
    tmp(:,2*j-1) = NumerDisp(2*k-1,n);
    tmp(:,2*j  ) = NumerDisp(2*k,  n);
    k = TriMesh.elem2edge(:,j);
    tmp(:,6+j) = NumerDisp(2*NumNds+k,n);
  end
  EltDsplDivOld = sum(squeeze(EltMatP0DivBR1).*tmp,2);
  GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) + alpha * EltDsplDivOld;

  % For solid/elas.: Incorporating Neumann (traction) boundary conditions 
  if NumNeumannEgsElas>0
    % For end-nodes 
    NeumannBC = zeros(NumNeumannEgsElas,2,2);
    j1 = TriMesh.edge(NeumannEdgeElas,1);  
    j2 = TriMesh.edge(NeumannEdgeElas,2);
    x1 = TriMesh.node(j1,1);  y1 = TriMesh.node(j1,2);
    x2 = TriMesh.node(j2,1);  y2 = TriMesh.node(j2,2);
    LenEgNeumann = sqrt((x2-x1).^2+(y2-y1).^2);
    NumQuadPts = size(GAUSSQUAD.LINE,1);
    for k=1:NumQuadPts
      qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
      tNval = EqnBIC.fxntN(qp,tnew);
      NeumannBC(:,1,1) = NeumannBC(:,1,1) + GAUSSQUAD.LINE(k,3)...
        * GAUSSQUAD.LINE(k,1) * tNval(:,1) .* LenEgNeumann;
      NeumannBC(:,1,2) = NeumannBC(:,1,2) + GAUSSQUAD.LINE(k,3)...
        * GAUSSQUAD.LINE(k,1) * tNval(:,2) .* LenEgNeumann;
      NeumannBC(:,2,1) = NeumannBC(:,2,1) + GAUSSQUAD.LINE(k,3)...
        * GAUSSQUAD.LINE(k,2) * tNval(:,1) .* LenEgNeumann;
      NeumannBC(:,2,2) = NeumannBC(:,2,2) + GAUSSQUAD.LINE(k,3)...
        * GAUSSQUAD.LINE(k,2) * tNval(:,2) .* LenEgNeumann;
    end
    % JL20171231: TO BE REVISED 
    for ig=1:NumNeumannEgsElas 
      j1 = TriMesh.edge(NeumannEdgeElas(ig),1);
      j2 = TriMesh.edge(NeumannEdgeElas(ig),2);
      GlbRhsDspl(2*j1-1) = GlbRhsDspl(2*j1-1) - NeumannBC(ig,1,1);
      GlbRhsDspl(2*j1  ) = GlbRhsDspl(2*j1  ) - NeumannBC(ig,1,2);
      GlbRhsDspl(2*j2-1) = GlbRhsDspl(2*j2-1) - NeumannBC(ig,2,1);
      GlbRhsDspl(2*j2  ) = GlbRhsDspl(2*j2  ) - NeumannBC(ig,2,2);
    end
    % For edges 
    NeumannBC = zeros(NumNeumannEgsElas,1);
    j1 = TriMesh.edge(NeumannEdgeElas,1);
    j2 = TriMesh.edge(NeumannEdgeElas,2);
    x1 = TriMesh.node(j1,1);  y1 = TriMesh.node(j1,2);
    x2 = TriMesh.node(j2,1);  y2 = TriMesh.node(j2,2);
    LenEgNeumann = sqrt((x2-x1).^2+(y2-y1).^2);
    nml = TriMesh.EgNml(NeumannEdgeElas,:);
    NumQuadPts = size(GAUSSQUAD.LINE,1);
    for k=1:NumQuadPts
      qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
      tNval = EqnBIC.fxntN(qp,tnew);
      psi = GAUSSQUAD.LINE(k,1) * GAUSSQUAD.LINE(k,2);
      NeumannBC = NeumannBC ...
        + GAUSSQUAD.LINE(k,3) * sum(tNval.*nml,2) * psi .* LenEgNeumann;
    end
    GlbRhsDspl(2*NumNds+NeumannEdgeElas) = ...
    GlbRhsDspl(2*NumNds+NeumannEdgeElas) - NeumannBC;
  end

  % For fluid/Darcy: Incorporating Neumann (flux) boundary conditions 
  if (NumNeumannEgsDarcy>0)
    NeumannBC = zeros(NumNeumannEgsDarcy,1);
    x1 = TriMesh.node(TriMesh.edge(NeumannEdgeDarcy,1),1);
    y1 = TriMesh.node(TriMesh.edge(NeumannEdgeDarcy,1),2);
    x2 = TriMesh.node(TriMesh.edge(NeumannEdgeDarcy,2),1);
    y2 = TriMesh.node(TriMesh.edge(NeumannEdgeDarcy,2),2);
    LenNeumannEg = sqrt((x2-x1).^2+(y2-y1).^2);
    NumQuadPts = size(GAUSSQUAD.LINE,1);
    for k=1:NumQuadPts
      qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
      uNval = EqnBIC.fxnuN(qp,tnew);
      NeumannBC = NeumannBC + GAUSSQUAD.LINE(k,3) * uNval .* LenNeumannEg; 
    end
    GlbRhsPres(NumEms+NeumannEdgeDarcy) = ...
    GlbRhsPres(NumEms+NeumannEdgeDarcy) - NeumannBC; 
  end

  % Assembling GlbRHS: Combination 
  GlbRHS = [GlbRhsDspl; GlbRhsPres];

  % For solution 
  sln = zeros(DOFsDspl+DOFsPres,1);

  
  % For Elas.: Enforcing Dirichlet boundary conditions as essential 
  % Assuming NumDirichletEgsElas>0 
  if (NumDirichletEgsElas>0) 
    % For nodes: 
    DirichletNdElas = TriMesh.node(DirichletNodeElas,:);
    uDval = EqnBIC.fxnuD(DirichletNdElas,tnew);
    sln(2*DirichletNodeElas-1) = uDval(:,1);
    sln(2*DirichletNodeElas  ) = uDval(:,2);
  end
  if (NumDirichletEgsElas>0)
    intgrl1 = zeros(NumDirichletEgsElas,1);
    intrpl = zeros(NumDirichletEgsElas,1);
    intgrl2 = zeros(NumDirichletEgsElas,1);
    x1 = TriMesh.node(TriMesh.edge(DirichletEdgeElas,1),1);
    y1 = TriMesh.node(TriMesh.edge(DirichletEdgeElas,1),2);
    x2 = TriMesh.node(TriMesh.edge(DirichletEdgeElas,2),1);
    y2 = TriMesh.node(TriMesh.edge(DirichletEdgeElas,2),2);
    LenDirichletEdgeElas = sqrt((x2-x1).^2+(y2-y1).^2);
    nml = TriMesh.EgNml(DirichletEdgeElas,:);
    NumQuadPts = size(GAUSSQUAD.LINE,1);
    for k=1:NumQuadPts
      qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
      uD = EqnBIC.fxnu(qp,tnew);
      psi = GAUSSQUAD.LINE(k,1)*GAUSSQUAD.LINE(k,2);
      intgrl1 = intgrl1 + GAUSSQUAD.LINE(k,3)*psi*LenDirichletEdgeElas;
      intrpl = GAUSSQUAD.LINE(k,1)*EqnBIC.fxnu(TriMesh.node(TriMesh.edge(DirichletEdgeElas,1),:),tnew)...
              + GAUSSQUAD.LINE(k,2)*EqnBIC.fxnu(TriMesh.node(TriMesh.edge(DirichletEdgeElas,2),:),tnew); 
      res = uD - intrpl;
      intgrl2 = intgrl2 + GAUSSQUAD.LINE(k,3)*sum(res.*nml,2).*LenDirichletEdgeElas;
    end
    sln(2*NumNds+DirichletEdgeElas) = intgrl2./intgrl1;
  end
  
    
  % For Darcy: Enforcing Dirichlet boundary conditions as essential
  % Assuming NumDirichletEgsDarcy>0
  if (NumDirichletEgsDarcy>0)
    pDavg = zeros(NumDirichletEgsDarcy,1);
    x1 = TriMesh.node(TriMesh.edge(DirichletEdgeDarcy,1),1);
    y1 = TriMesh.node(TriMesh.edge(DirichletEdgeDarcy,1),2);
    x2 = TriMesh.node(TriMesh.edge(DirichletEdgeDarcy,2),1);
    y2 = TriMesh.node(TriMesh.edge(DirichletEdgeDarcy,2),2);
    NumQuadPts = size(GAUSSQUAD.LINE,1);
    for k=1:NumQuadPts
      qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
      pDval = EqnBIC.fxnpD(qp,tnew);
      pDavg = pDavg + GAUSSQUAD.LINE(k,3) * pDval;
    end
    sln(DOFsDspl+NumEms+DirichletEdgeDarcy) = pDavg;
  end
  
  % Reducing...
  GlbRHS = GlbRHS - GlbMat*sln;
%   GlbRHS


  % For flag 
  flag = ones(DOFsDspl+DOFsPres,1);
  % Elas D.B.C.nodes
  flag(2*DirichletNodeElas-1) = 0;
  flag(2*DirichletNodeElas  ) = 0;
  % Elas D.B.C.edges
  flag(2*NumNds+DirichletEdgeElas) = 0;
  % Darcy D.B.C edges
  flag(DOFsDspl+NumEms+DirichletEdgeDarcy) = 0;
  FreeNdEmFreeEg = find(flag);
  
  % Solving the reduced linear system 
  sln(FreeNdEmFreeEg) = GlbMat(FreeNdEmFreeEg,FreeNdEmFreeEg) \ GlbRHS(FreeNdEmFreeEg);

  % Extracting and storing 
  NumerDisp(:,n+1) = sln(1:DOFsDspl);
  NumerPres(:,n+1) = sln(DOFsDspl+(1:DOFsPres));
  
  %% Auxiliary quantities 
    % Elementwise 4 coeffs. in WG bas.fxns. 
    cof = zeros(NumEms,4);
    cof(:,1) = sln(DOFsDspl+1:DOFsDspl+NumEms);
    for j=1:3
     cof(:,1+j) = sln(DOFsDspl+NumEms+TriMesh.elem2edge(:,j));
    end
    % Elementwise cofs. in RT0 bas.: Disc.wk.grad. of WG0 bas.fxns.
    % CDWGB = zeros(NumEms,4,3) 
    CDWGB = WG_TriP0P0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(TriMesh);
    CDWGS = zeros(NumEms,3);
    for j=1:3
     CDWGS(:,j) = sum(cof.*squeeze(CDWGB(:,:,j)),2);
    end
    % GM = zeros(NumEms,3,3);
    % GMK = zeros(NumEms,3,3);
    GM = Hdiv_TriRT0_NmlzBas_GramMat(TriMesh); 
    GMK = Hdiv_TriRT0_NmlzBas_GramMatK(TriMesh,PermK);
    % Projection coeffs. in RT0 nmlz. bas. for permeability * bas.fxn.
    ProjCof = zeros(NumEms,3,3);  % j-th fxn. projected onto k-th fxn. 
    for ie=1:NumEms  % TO BE REVISED FOR EFFECIENCY! 
     EltGM = squeeze(GM(ie,:,:));
     for j=1:3
      EltRHS = squeeze(GMK(ie,:,j));
      ProjCof(ie,j,:) = EltGM\EltRHS';
     end
    end

    %% Numerical velocity: Coeffs. in RT elementwise 
    for k=1:3
    NumerVelCof(:,k,n) = -sum(CDWGS.*squeeze(ProjCof(:,:,k)),2);
    end
end
fprintf('\n');

return;