%% run_StokesDarcy_QuadriCGBR1WGAC0_3domH.m
% Stokes-Darcy coupling: 3 domains horizontally: Stokes/Darcy/Stokes 
% -- Stokes solved by CG(BR1,Q0); 
% -- Darcy solved by WG(P0,P0;AC0);
% James Liu, ColoState; Graham Harper, Sandia; 2019/07--2021/04 

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

%% Setting a problem: Domain1/2/3: Stokes/Darcy/Stokes, horizontal lineup 
% Boundary condition types: C, S,D,S, S, S,D,S, S (totaly 8 sides) 
% 11 for interface#1, 12 for interface#2 
% Check to make sure: xb1=xa2, xb2=xa3, yc1=yc2=yc3, yd1=yd2=yd3 
% 
% EqnBC1 = EqnBC_StokesDarcy3domH_LiuEx_Part1Stokes; 
% EqnBC2 = EqnBC_StokesDarcy3domH_LiuEx_Part2Darcy; 
% EqnBC3 = EqnBC_StokesDarcy3domH_LiuEx_Part3Stokes; 
% BndryCondType = [0; 1;1;1; 1; 1;1;1; 1];
% xa1 = -1;  xb1 = 0;  yc1 = 0;  yd1 = 1;
% xa2 =  0;  xb2 = 1;  yc2 = 0;  yd2 = 1;
% xa3 =  1;  xb3 = 2;  yc3 = 0;  yd3 = 1;
% 
EqnBC1 = EqnBC_StokesDarcy3domH_GrahamExVar_Part1Stokes; 
EqnBC2 = EqnBC_StokesDarcy3domH_GrahamExVar_Part2Darcy; 
EqnBC3 = EqnBC_StokesDarcy3domH_GrahamExVar_Part3Stokes; 
BndryCondType = [0; 1;2;1; 2; 1;2;1; 1];
xa1 = -1;  xb1 = 0;  yc1 = 0;  yd1 = 1;
xa2 =  0;  xb2 = 1;  yc2 = 0;  yd2 = 1;
xa3 =  1;  xb3 = 2;  yc3 = 0;  yd3 = 1;
% 
BndryCondType1 = [0; BndryCondType(2); 11; BndryCondType(8); BndryCondType(9)];
BndryCondType2 = [0; BndryCondType(3); 12; BndryCondType(7); 11];
BndryCondType3 = [0; BndryCondType(4); BndryCondType(5); BndryCondType(6); 12];
BndryDescMat1 = [xa1, yc1, xb1, yc1,  0,-1; ... 
                 xb1, yc1, xb1, yd1,  1, 0; ... 
                 xb1, yd1, xa1, yd1,  0, 1; ... 
                 xa1, yd1, xa1, yc1, -1, 0]; 
BndryDescMat2 = [xa2, yc2, xb2, yc2,  0,-1; ... 
                 xb2, yc2, xb2, yd2,  1, 0; ... 
                 xb2, yd2, xa2, yd2,  0, 1; ... 
                 xa2, yd2, xa2, yc2, -1, 0]; 
BndryDescMat3 = [xa3, yc3, xb3, yc3,  0,-1; ... 
                 xb3, yc3, xb3, yd3,  1, 0; ... 
                 xb3, yd3, xa3, yd3,  0, 1; ... 
                 xa3, yd3, xa3, yc3, -1, 0];

%% Preparing meshes 
% Generating three indepdent rectangular meshes, need check ny1=ny2=ny3 
n = 16;  delta = 0;  status = 2;
nx1 = n;  ny1 = n;  
QuadriMesh1 = RectDom_QuadriMesh_GenLogiRectPtbSmth(xa1,xb1,nx1, yc1,yd1,ny1, delta, status);
nx2 = n;  ny2 = n;  
QuadriMesh2 = RectDom_QuadriMesh_GenLogiRectPtbSmth(xa2,xb2,nx2, yc2,yd2,ny2, delta, status);
nx3 = n;  ny3 = n;  
QuadriMesh3 = RectDom_QuadriMesh_GenLogiRectPtbSmth(xa3,xb3,nx3, yc3,yd3,ny3, delta, status);
% 
QuadriMesh1 = QuadriMesh_Enrich2(QuadriMesh1, BndryDescMat1);
QuadriMesh1 = QuadriMesh_Enrich3(QuadriMesh1, BndryDescMat1);
QuadriMesh1 = QuadriMesh_Enrich4(QuadriMesh1, BndryDescMat1);
% 
QuadriMesh2 = QuadriMesh_Enrich2(QuadriMesh2, BndryDescMat2);
QuadriMesh2 = QuadriMesh_Enrich3(QuadriMesh2, BndryDescMat2);
QuadriMesh2 = QuadriMesh_Enrich4(QuadriMesh2, BndryDescMat2);
% 
QuadriMesh3 = QuadriMesh_Enrich2(QuadriMesh3, BndryDescMat3);
QuadriMesh3 = QuadriMesh_Enrich3(QuadriMesh3, BndryDescMat3);
QuadriMesh3 = QuadriMesh_Enrich4(QuadriMesh3, BndryDescMat3);
% Concatenation 
QuadriMesh = QuadriMesh_ConcatHoriz31(QuadriMesh1,QuadriMesh2,QuadriMesh3);
% if QuadriMesh2.NumEms<=2^8; show_QuadriMesh_labels(QuadriMesh1,12,1,1,1); end
% if QuadriMesh1.NumEms<=2^8; show_QuadriMesh_labels(QuadriMesh1,11,1,1,1); end
% if QuadriMesh3.NumEms<=2^8; show_QuadriMesh_labels(QuadriMesh1,13,1,1,1); end
% if QuadriMesh.NumEms<=2^8; show_QuadriMesh_labels(QuadriMesh,14,1,0,1); end

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

%% Assemling and Solving 
[sln] = StokesDarcy_QuadriCGBR1WGAC0_3domH_AsmSlv( ... 
  EqnBC1, BndryCondType1, QuadriMesh1, ... 
  EqnBC2, BndryCondType2, QuadriMesh2, ... 
  EqnBC3, BndryCondType3, QuadriMesh3, GAUSSQUAD);

%% Postprocessing for presentation 
[NumerPresEmAvg, NumerVelEmAvg] ... 
  = StokesDarcy_QuadriCGBR1WGAC0_3domH_Presn( ... 
    sln, QuadriMesh1, QuadriMesh2, QuadriMesh3, GAUSSQUAD, EqnBC2);

%% Presenting numerical and graphical results 
% The old way 
% show_QuadriMesh_ScaVecEm_mix(QuadriMesh, NumerPresEmAvg, NumerVelEmAvg, ...
%   27, 'Stokes-Darcy-Stokes: Numerical pressure & velocity', 1);
% 
DOFs1Stokes = 2*QuadriMesh1.NumNds + QuadriMesh1.NumEgs + QuadriMesh1.NumEms;
DOFs3Stokes = 2*QuadriMesh3.NumNds + QuadriMesh3.NumEgs + QuadriMesh3.NumEms;
DOFs2Darcy  = QuadriMesh2.NumEms + QuadriMesh2.NumEgs;
slnDarcy = sln(DOFs1Stokes+(1:DOFs2Darcy));
slnDarcyEm = slnDarcy(1:QuadriMesh2.NumEms);
reshape(slnDarcyEm,n,n);
slnDarcyEg = slnDarcy(QuadriMesh2.NumEms+(1:QuadriMesh2.NumEgs));
slnDarcyEg(1:n);
PermK = Darcy_SmplnPerm_QuadriMesh(EqnBC2.fxnK, QuadriMesh2, GAUSSQUAD);
CofQhg = zeros(QuadriMesh2.NumEms,4);
% 
[NumerPresEmAvg2, NumerVelEmAvg2] = ... 
  DarcyGen_WG_QuadriP0P0AC0_PresVelFlux( ...
    EqnBC2, QuadriMesh2, PermK, slnDarcy, CofQhg, GAUSSQUAD);
% 
% [NumerPresEmAvg2, NumerVelEmAvg2, NumerVelNd2] = ... 
%   DarcyGen_WG_QuadriP0P0AC0_Presn( ...
%     EqnBC2, QuadriMesh2, PermK, slnDarcy, CofQhg, GAUSSQUAD);
% 
% [NumerPresEm,NumerVelCofRT0,~,~,~] = ...
%   Darcy_WG_QuadriQ0Q0RT0_PresVelFlux(QuadriMesh2,PermK,slnDarcy,GAUSSQUAD);
% NumerPresEmAvg2 = NumerPresEm;
% NumerVelEmAvg2 = NumerVelCofRT0(:,[1,2]);
% 
show_QuadriMesh_ScaVecEm_mix(QuadriMesh2, NumerPresEmAvg2, NumerVelEmAvg2, ...
  903, 'Darcy: Numerical pressure & velocity', 1);
% 
NumerVelNd1 = [sln(1:2:(2*QuadriMesh1.NumNds)) sln(2:2:(2*QuadriMesh1.NumNds))];
NumerPresEmAvg1 = NumerPresEmAvg(1:QuadriMesh1.NumEms);
% 
% NumerPresEmAvg2 = NumerPresEmAvg(QuadriMesh1.NumEms + (1:QuadriMesh2.NumEms));
% NumerVelEmAvg2 = NumerVelEmAvg(QuadriMesh1.NumEms +   (1:QuadriMesh2.NumEms),:);
% 
NumerVelNd3 = [sln(DOFs1Stokes + DOFs2Darcy + (1:2:(2*QuadriMesh3.NumNds))) ...
               sln(DOFs1Stokes + DOFs2Darcy + (2:2:(2*QuadriMesh3.NumNds)))];
NumerPresEmAvg3 = NumerPresEmAvg(QuadriMesh1.NumEms + QuadriMesh2.NumEms + (1:QuadriMesh3.NumEms));
% 
% plot all on figure 27 
% (it must be this order or velocity will not display properly)
show_QuadriMesh_ScaEmVecNd_mix(QuadriMesh3, NumerPresEmAvg3, NumerVelNd3, ...
  27, 'Stokes-Darcy-Stokes: Numerical pressure & velocity', 1);
show_QuadriMesh_ScaVecEm_mix(QuadriMesh2, NumerPresEmAvg2, NumerVelEmAvg2, ...
  27, 'Stokes-Darcy-Stokes: Numerical pressure & velocity', 1);
show_QuadriMesh_ScaEmVecNd_mix(QuadriMesh1, NumerPresEmAvg1, NumerVelNd1, ...
  27, 'Stokes-Darcy-Stokes: Numerical pressure & velocity', 1);
% 
plot([0,0],[0,1], 'k-')  % plot left interface
plot([1,1],[0,1], 'k-')  % plot right interface
axis tight;
% 
ax = gca; % grab figure axes
ax.TitleFontSizeMultiplier = 2; % double title size so I can take a large screenshot
% 
colBar = colorbar();
colBar.Limits(1) = 0; % set colorbar limits to start at 0

% %% For debugging 
% ExactPresEmAvg1 = EqnBC1.fxnp(QuadriMesh1.EmCntr);
% ExactPresEmAvg2 = EqnBC2.fxnp(QuadriMesh2.EmCntr);
% ExactPresEmAvg3 = EqnBC3.fxnp(QuadriMesh3.EmCntr);
% ExactVelEmAvg1 = EqnBC1.fxnu(QuadriMesh1.EmCntr);
% ExactVelEmAvg2 = EqnBC2.fxnu(QuadriMesh2.EmCntr);
% ExactVelEmAvg3 = EqnBC3.fxnu(QuadriMesh3.EmCntr);
% ExactPresEmAvg = [ExactPresEmAvg1; ExactPresEmAvg2; ExactPresEmAvg3];
% ExactVelEmAvg = [ExactVelEmAvg1; ExactVelEmAvg2; ExactVelEmAvg3];
% show_QuadriMesh_ScaVecEm_mix(QuadriMesh, ExactPresEmAvg, ExactVelEmAvg, ...
%   93, 'Stokes-Darcy-Stokes: Exact pressure & velocity', 1);
% % max(max(NumerVelEmAvg-ExactVelEmAvg))
% % DOFs1Stokes = 2*QuadriMesh1.NumNds + QuadriMesh1.NumEgs + QuadriMesh1.NumEms;
% % DOFs3Stokes = 2*QuadriMesh3.NumNds + QuadriMesh3.NumEgs + QuadriMesh3.NumEms;
% % DOFs2Darcy  = QuadriMesh2.NumEms + QuadriMesh2.NumEgs;
% % DOFs = DOFs1Stokes + DOFs2Darcy + DOFs3Stokes;
% % [DOFs1Stokes,DOFs2Darcy,DOFs3Stokes,DOFs]
% % figure(97);  spy(GlbMat);
% % figure(98);  plot(GlbRHS);
% % figure(99);  plot(sln);

return;