%% run_StokesDarcy_QuadriCGBR1WGAC0_1meshBT.m 
% The 1-(quadri.)mesh approach with Darcy at bottom and Stokes at top 
% -- Stokes solved by CG(BR1,Q0); 
% -- Darcy solved by WG(P0,P0;AC0);
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% Graham Harper, James Liu, ColoState; 2019/07--2021/04 

% clc;  
clear all;  close all;
format compact;  format short e;

%% Setting up a problem to be solved 
% Boundary condition type: 1 Dirichlet, 2 Neumann, 0 convenience 
% A domain & its boundary (start/ending points, outward unit normal vec.) 
% Stokes/Darcy domains split in the middle 
% 
% EqnBC = EqnBC_StokesDarcy_ChenWangWangIMAJNA2016_Ex; 
% BndryCondType = [0; 1; 1; 1; 1; 1; 1];  % [C; D; D; S; S; S; D]
% xa = 0;  xb = pi;  yc = -1;  yd = 1; 
% 
% EqnBC = EqnBC_StokesDarcy_DiscaQuarteRevMatComplut2009_ExP383Var;
% BndryCondType = [0; 1; 2; 1; 2; 1; 2];  % [C; D; D; S; S; S; D]
% xa = 0;  xb = 1;  yc = 0;  yd = 2; 
% 
EqnBC = EqnBC_StokesDarcy_Harper3CMAME2021_Ex3;
BndryCondType = [0; 2;2; 1;1;1; 2];  % [C; D; D; S; S; S; D]
xa = 0;  xb = 2;  yc = -1;  yd = 1;
% 
BndryDescMat = [xa, yc,        xb, yc,         0,-1; ... 
                xb, yc,        xb, (yc+yd)/2,  1, 0; ... 
                xb, (yc+yd)/2, xb, yd,         1, 0; ... 
                xb, yd,        xa, yd,         0, 1; ... 
                xa, yd,        xa, (yc+yd)/2, -1, 0; ...
                xa, (yc+yd)/2, xa, yc,        -1, 0];

%% Discretization: Mesh generation and preparation for the 1-domain 
% For a rectangular or smooth quadrilateral mesh 
% For EqnBC_StokesDarcy_Harper3CMAME2021_Ex3, n = multiple of 20 
n = 40;  nx = n;  ny= n;  delta = 0;  status = 2;
QuadriMesh = RectDom_QuadriMesh_GenLogiRectPtbSmth(...
  xa,xb,nx,yc,yd,ny,delta,status);
% % For a trapezoidal mesh proposed in ArnoldBoffiFalk_MCOM_2002 
% n = 2^4;  nx = n;  ny= n;  slant = 0.35;
% QuadriMesh = RectDom_QuadriMesh_GenTrapezoid(xa,xb,nx,yc,yd,ny,slant,2);
QuadriMesh = QuadriMesh_Enrich2(QuadriMesh,BndryDescMat);
QuadriMesh = QuadriMesh_Enrich3(QuadriMesh,BndryDescMat);
QuadriMesh = QuadriMesh_Enrich4(QuadriMesh,BndryDescMat);
if QuadriMesh.NumEms<=2^8
  show_QuadriMesh_labels(QuadriMesh,11,0,0,0);
end
% Enriching mesh info for the (1-domain) Stokes-Darcy (top-bottom) coupling
QuadriMesh = QuadriMesh_Enrich_StokesDarcy(QuadriMesh);

%% Setting up Gaussian quadratures 
GAUSSQUAD = SetGaussQuad(5,25,13);

%% Stokes-Darcy solved by CG(BR1,Q0) for Stokes and WG(P0,P0;AC0) for Darcy 
[sln,CofQhg] = StokesDarcy_QuadriCGBR1WGAC0_1meshBT_AsmSlv( ... 
  EqnBC, QuadriMesh, BndryCondType, GAUSSQUAD);

%% Computing errors if exact solutions exist 
[hErrStokesDarcy,L2ErrVelStokes, L2ErrPresDarcy,L2ErrVelDarcy,L2ErrDivDarcy] ... 
  = StokesDarcy_QuadriCGBR1WGAC0_1meshBT_Err( ...
    EqnBC, QuadriMesh, sln, GAUSSQUAD);

%% Postprocessing for presentation 
[NumerVelEmAvg, NumerPresEmAvg] = ... 
  StokesDarcy_QuadriCGBR1WGAC0_1meshBT_Presn( ...
  EqnBC, QuadriMesh, sln, CofQhg, GAUSSQUAD);

%% Presenting numerical and graphical results 
% shift = sum(NumerPresEmAvg.*QuadriMesh.area)/sum(QuadriMesh.area);
% NumerPresEmAvg = NumerPresEmAvg - shift;
% [hErrStokesDarcy,L2ErrVelStokes,L2ErrPresDarcy]
[L2ErrVelStokes,L2ErrPresDarcy,L2ErrVelDarcy,L2ErrDivDarcy]
%
ExactPresEmCntr = zeros(QuadriMesh.NumEms,1);
ExactPresEmCntr(QuadriMesh.StokesEms) = ... 
  EqnBC.fxnp_Stokes(QuadriMesh.EmCntr(QuadriMesh.StokesEms,:));
ExactPresEmCntr(QuadriMesh.DarcyEms) = ... 
  EqnBC.fxnp_Darcy(QuadriMesh.EmCntr(QuadriMesh.DarcyEms,:));
% 
ExactVelEmCntr = zeros(QuadriMesh.NumEms,2);
ExactVelEmCntr(QuadriMesh.StokesEms,:) = ... 
  EqnBC.fxnu_Stokes(QuadriMesh.EmCntr(QuadriMesh.StokesEms,:));
ExactVelEmCntr(QuadriMesh.DarcyEms,:) = ... 
  EqnBC.fxnu_Darcy(QuadriMesh.EmCntr(QuadriMesh.DarcyEms,:));
show_QuadriMesh_ScaVecEm_mix(QuadriMesh, ExactPresEmCntr, ExactVelEmCntr, ...
  23, 'Stokes-Darcy: Exact pressure & velocity', 1);
figure(23);  axis tight;  axis square;
show_QuadriMesh_ScaVecEm_mix(QuadriMesh, NumerPresEmAvg, NumerVelEmAvg, ...
  27, 'Stokes-Darcy: Numerical pressure & velocity', 2);
figure(27);  axis tight;  axis square;
% print -f23 -dpng Ex_fig_StokesDarcy_ExactVelPres_n.png
% print -f27 -dpng Ex_fig_StokesDarcy_NumerVelPres_n.png

return;