%% JL20200711: BEING REVISED 
function [hErrStokesDarcy,L2ErrVelStokes, ...
  L2ErrPresDarcy,L2ErrVelDarcy,L2ErrDivDarcy] ... 
  = StokesDarcy_QuadriCGBR1WGAC0_1meshBT_Err( ...
    EqnBC, QuadriMesh, sln, GAUSSQUAD) 
% The 1-(quadri.)mesh approach with Darcy at bottom and Stokes at top 
% -- QuadriMesh is 1-mesh for the concatenation of two domains; 
% -- sln contains all lin.comb.coeffs. for bas.fxns. Stokes & Darcy 
% Graham Harper, James Liu, ColoState; 2019/07--2020/07 

%% Mesh info 
% NumNds = QuadriMesh.NumNds;
NumEgs = QuadriMesh.NumEgs;
NumEms = QuadriMesh.NumEms;
% area = QuadriMesh.area;
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% More 
LblEmDarcy  = QuadriMesh.DarcyEms;
LblEmStokes = QuadriMesh.StokesEms;
NumEmsDarcy = max(size(QuadriMesh.DarcyEms));
NumEgsDarcy = nnz(QuadriMesh.posDarcyEg);
NumDOFsDarcy = NumEmsDarcy + NumEgsDarcy;
NumNdsStokes = nnz(QuadriMesh.posStokesNd);
NumEgsStokes = nnz(QuadriMesh.posStokesEg);
NumEmsStokes = length(LblEmStokes);
NumDOFsStokes = 2*NumNdsStokes + NumEgsStokes + NumEmsStokes;
DOFs = NumDOFsDarcy + NumDOFsStokes;

%% JL20191107: UNDER THE ASSUMPTION "Darcy elements go first" 
%% Sampling the permeability (2x2 SPD matrix)
PermK = Darcy_SmplnPerm_QuadriMesh(EqnBC.fxnK_Darcy, QuadriMesh, GAUSSQUAD);

%% Stokes: L2-norm of errors in velocity 
% Note the tricks for retrieving cof. of Stokes numerical velocity 
cof = zeros(NumEms,12);
for j=1:4  % Treating vertices and edges in one loop 
  k = QuadriMesh.posStokesNd(QuadriMesh.elem(QuadriMesh.StokesEms,j));
  cof(QuadriMesh.StokesEms,2*j-1) = sln(NumDOFsDarcy+2*k-1);
  cof(QuadriMesh.StokesEms,2*j  ) = sln(NumDOFsDarcy+2*k  );
  k = QuadriMesh.posStokesEg(QuadriMesh.elem2edge(QuadriMesh.StokesEms,j));
  cof(QuadriMesh.StokesEms,8+j) = ... 
    sln(NumDOFsDarcy+2*NumNdsStokes+k) .* QuadriMesh.EmEgSign(QuadriMesh.StokesEms,j);
end
% Computing errors in L2-norm 
EmErrVel = zeros(NumEms,1);
NumerVelEm = zeros(NumEms,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
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);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  u = EqnBC.fxnu_Stokes(qp);
  BFVal = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat);
  NumerVelEm(:,1) = sum(cof.*squeeze(BFVal(:,:,1)),2);
  NumerVelEm(:,2) = sum(cof.*squeeze(BFVal(:,:,2)),2);
  uerr = u - NumerVelEm;
  EmErrVel = EmErrVel + (uerr(:,1).^2+uerr(:,2).^2)...
    .* jac * GAUSSQUAD.RECT(k,3);
end
L2ErrVelStokes = sqrt(sum(EmErrVel(QuadriMesh.StokesEms)));

%% Darcy: L2-norm of errors in pressure 
% Note the tricks here for retrieving Darcy numerical pressure 
NumerPresEm = zeros(NumEms,1);
NumerPresEm(QuadriMesh.DarcyEms) = sln(1:NumEmsDarcy);
% Computing errors in L2-norm 
EmErrPres = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  pval = EqnBC.fxnp_Darcy(qp);
  perr = pval - NumerPresEm;
  EmErrPres = EmErrPres + (perr.^2) .* jac * GAUSSQUAD.RECT(k,3);
end
L2ErrPresDarcy = sqrt(sum(EmErrPres(QuadriMesh.DarcyEms)));

%% Darcy: Auxialiary 
% Note the tricks here for retrieving cofs. for Darcy  
cof = zeros(NumEms,5);  % Elementwise 5 coeffs. in WG bas.fxns. 
cof(QuadriMesh.DarcyEms,1) = sln(1:NumEmsDarcy);
for j=1:4 
  k= QuadriMesh.posDarcyEg(QuadriMesh.elem2edge(QuadriMesh.DarcyEms,j));
  cof(QuadriMesh.DarcyEms,1+j) = sln(NumEmsDarcy+k);
end
% Elementwise cofs. in AC0 nmlz.Piola.bas.: Disc.wk.grad. of WG0 bas.fxns.
% CDWGB = zeros(NumEms,5,4) 
CDWGB = WG_QuadriP0P0AC0_CofAC0NmlzPiolaBas_DiscWkGradBasFxn(QuadriMesh,GAUSSQUAD);
CDWGS = zeros(NumEms,4);
for j=1:4
  CDWGS(:,j) = sum(cof.*squeeze(CDWGB(:,:,j)),2);
end
% GM = zeros(NumEms,4,4);
% GMK = zeros(NumEms,4,4);
GM = Hdiv_QuadriAC0_NmlzPiolaBas_GramMat(QuadriMesh,GAUSSQUAD);
GMK = Hdiv_QuadriAC0_NmlzPiolaBas_GramMatK(QuadriMesh,PermK,GAUSSQUAD);
% Projection coeffs. in AC0 nmlz.Piola.bas. for permeability * bas.fxn.
ProjCof = zeros(NumEms,4,4);  % j-th fxn. projected onto k-th fxn. 
for ie=1:NumEms  % TO BE REVISED FOR EFFECIENCY! 
  EltGM = squeeze(GM(ie,:,:));
  for j=1:4
    EltRHS = squeeze(GMK(ie,:,j));
    ProjCof(ie,j,:) = EltGM\EltRHS';
  end
end
% Numerical velocity: Coeffs. in AC0 elementwise 
NumerVelCofAC0 = zeros(NumEms,4);
for k=1:4
  NumerVelCofAC0(:,k) = -sum(CDWGS.*squeeze(ProjCof(:,:,k)),2);
end

%% Darcy: L2-norm of errors in velocity 
% Using NumerVelCofAC0 
MatJ = zeros(NumEms,2,2);
MatP = zeros(NumEms,2,2);
ErrVelEm = zeros(NumEms,1);
NV = zeros(NumEms,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  J11 = CofA(:,2) + CofA(:,4)*yhat;  J12 = CofA(:,3) + CofA(:,4)*xhat;
  J21 = CofB(:,2) + CofB(:,4)*yhat;  J22 = CofB(:,3) + CofB(:,4)*xhat;
  jac = J11.*J22 - J12.*J21;
  MatJ(:,1,1) = J11;  MatJ(:,1,2) = J12;
  MatJ(:,2,1) = J21;  MatJ(:,2,2) = J22;
  for i=1:2
    for j=1:2
      MatP(:,i,j) = MatJ(:,i,j)./jac;
    end
  end
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  X = qp(:,1) - QuadriMesh.EmCntr(:,1);
  Y = qp(:,2) - QuadriMesh.EmCntr(:,2);
  U = MatP(:,1,1) .* xhat - MatP(:,1,2) .* yhat;
  V = MatP(:,2,1) .* xhat - MatP(:,2,2) .* yhat;
  NV(:,1) = NumerVelCofAC0(:,1) + NumerVelCofAC0(:,3).*X + NumerVelCofAC0(:,4).*U;
  NV(:,2) = NumerVelCofAC0(:,2) + NumerVelCofAC0(:,3).*Y + NumerVelCofAC0(:,4).*V;
  uval = EqnBC.fxnu_Darcy(qp);
  uerr = uval - NV;
  ErrVelEm = ErrVelEm + (uerr(:,1).^2+uerr(:,2).^2) .*jac * GAUSSQUAD.RECT(k,3);
end
L2ErrVelDarcy = sqrt(sum(ErrVelEm(QuadriMesh.DarcyEms)));

%% Darcy: L2-norm of errors in velocity div  
% Note: Exact divergence = source: div(u) = f in Darcy equation 
ErrDivEm = zeros(NumEms,1);
NumerDivEm = 2*NumerVelCofAC0(:,3);  % Numerical velocity divergence 
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat) ...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat) ;
  divu = EqnBC.fxns_Darcy(qp);
  derr = divu - NumerDivEm;
  ErrDivEm = ErrDivEm + (derr.^2) .* jac * GAUSSQUAD.RECT(k,3);
end
L2ErrDivDarcy = sqrt(sum(ErrDivEm(QuadriMesh.DarcyEms)));

%% JL20200303: TO BE FINISHED BY Hraham 
%% Computing Stokes interpolation coefficients 

%% Computing Darcy projection to P0 
% For elements 
DarcyPresProjP0Em = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  pval = EqnBC.fxnp_Darcy(qp);
  DarcyPresProjP0Em = DarcyPresProjP0Em + pval.*jac*GAUSSQUAD.RECT(k,3);
end
% For edges 
DarcyPresProjP0Eg = zeros(NumEgs,1);
x1 = QuadriMesh.node(QuadriMesh.edge(:,1),1);
y1 = QuadriMesh.node(QuadriMesh.edge(:,1),2);
x2 = QuadriMesh.node(QuadriMesh.edge(:,2),1);
y2 = QuadriMesh.node(QuadriMesh.edge(:,2),2);
NumQuadPts = size(GAUSSQUAD.LINE,1);
for k=1:NumQuadPts
  qp = [GAUSSQUAD.LINE(k,1)*x1 + GAUSSQUAD.LINE(k,2)*x2,...
        GAUSSQUAD.LINE(k,1)*y1 + GAUSSQUAD.LINE(k,2)*y2];
  pval = EqnBC.fxnp_Darcy(qp);  
  DarcyPresProjP0Eg = DarcyPresProjP0Eg + pval * GAUSSQUAD.LINE(k,3);
end

%% Computing h-norm for Darcy discrete error 
CofDarcyPresEm = zeros(NumEms,5);
CofDarcyPresEm(LblEmDarcy,1) = ...
  DarcyPresProjP0Em(LblEmDarcy) - sln(LblEmDarcy);
for j=1:4
  k = QuadriMesh.posDarcyEg(QuadriMesh.elem2edge(LblEmDarcy,j));
  CofDarcyPresEm(LblEmDarcy,1+j) = ... 
    DarcyPresProjP0Eg(QuadriMesh.elem2edge(LblEmDarcy,j)) - sln(NumEmsDarcy+k);
end
GMK = Hdiv_QuadriAC0_NmlzPiolaBas_GramMatK(QuadriMesh, PermK, GAUSSQUAD);
CDWGB = WG_QuadriP0P0AC0_CofAC0NmlzPiolaBas_DiscWkGradBasFxn(... 
  QuadriMesh, GAUSSQUAD);
cof = zeros(NumEms,4);
for k=1:4
  for j=1:5
    cof(:,k) = cof(:,k) + CofDarcyPresEm(:,j).*CDWGB(:,j,k);
  end
end
hErrDarcyEm = zeros(NumEms,1);
for kr=1:4
  for kc=1:4
    hErrDarcyEm = hErrDarcyEm + cof(:,kr).*GMK(:,kr,kc).*cof(:,kc);
  end
end
hErrDarcy = sqrt(sum(hErrDarcyEm(LblEmDarcy)));

% %% Computing h-norm for Stokes-Darcy combined  
% hErrStokesDarcy = sqrt(hErrStokes^2+hErrDarcy^2);
hErrStokesDarcy = sqrt(hErrDarcy^2);

return;