function [NumerPresEmAvg, NumerVelEmAvg] ... 
  = StokesDarcy_QuadriCGBR1WGAC0_3domH_Presn( ... 
    sln, QuadriMesh1, QuadriMesh2, QuadriMesh3, GAUSSQUAD, EqnBC2) 
%% Presentation 
% James Liu, ColoState; 2019/07--2020/09 

%% Auxiliary stuff 
DOFs1Stokes = 2*QuadriMesh1.NumNds + QuadriMesh1.NumEgs + QuadriMesh1.NumEms;
DOFs3Stokes = 2*QuadriMesh3.NumNds + QuadriMesh3.NumEgs + QuadriMesh3.NumEms;
DOFs2Darcy  = QuadriMesh2.NumEms + QuadriMesh2.NumEgs;

%% Extracting solutions for subdomains 
sln1 = sln([1:DOFs1Stokes]);
sln2 = sln(DOFs1Stokes+[1:DOFs2Darcy]);
sln3 = sln(DOFs1Stokes+DOFs2Darcy+[1:DOFs3Stokes]);

%% JL20200815: TO BE REVISED FOR Domain #2 (Darcy) 
%% Computing presentation quantities for subdomains 
[NumerVelEmAvg1, NumerDivEmAvg1, NumerPresEmAvg1, NumerVelNd1] ... 
  = Stokes_CG_QuadriBR1Q0_Presn(QuadriMesh1, sln1, GAUSSQUAD);
[NumerVelEmAvg3, NumerDivEmAvg3, NumerPresEmAvg3, NumerVelNd3] ... 
  = Stokes_CG_QuadriBR1Q0_Presn(QuadriMesh3, sln3, GAUSSQUAD);
% Awkward old code for Darcy 
PermK = Darcy_SmplnPerm_QuadriMesh(EqnBC2.fxnK, QuadriMesh2, GAUSSQUAD); 
CofQhg = zeros(QuadriMesh1.NumEms,4);
[NumerPresCofP0,NumerVelCofAC0] = ... 
  DarcyGen_WG_QuadriP0P0AC0_PresVelFlux( ...
    EqnBC2, QuadriMesh2, PermK, sln2, CofQhg, GAUSSQUAD);
NumerPresEmAvg2 = NumerPresCofP0(1:QuadriMesh2.NumEms);
NumerVelEmCntr2 = NumerVelCofAC0(:,1:2);

%% Now merging 
NumerPresEmAvg = [NumerPresEmAvg1; NumerPresEmAvg2; NumerPresEmAvg3];
NumerVelEmAvg = [NumerVelEmAvg1; NumerVelEmCntr2; NumerVelEmAvg3];

return;