%% run_LinPoroElas_MS_QuadriIE_WGP02P02AC02P0_WGP0P0AC0.m 
% Solving linear poroelasticity as a monolithic system (MS) 
% on a quadrilateral mesh with implicit Euler (IE) and 
% -- WG(P0^2,P0^2;AC0^2,P0) for lin.elasticity (solid displacement)
% -- WG(P0,P0;AC0) for Darcy flow (fluid pressure) 
% -- No need for penalty 
% First elas. then Darcy, since pressure changes faster 
% James Liu, ColoState; Zhuoran Wang, SYSU; 2020/09--2021/04 

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

%% Setting up the problem to be solved 
% Boundary condition types: 0 convenience, 1 Dirichlet, 2 Neumann 
% EqnBIC = EqnBIC_LinPoroElas_KansRivieJSC2018_ExVar;
EqnBIC = EqnBIC_LinPoroElas_LiuWangICCS2021_Ex1;
BndryCondTypeElas  = [0;1;1;1;1];
BndryCondTypeDarcy = [0;1;1;1;1];
xa = 0;  xb = 1;  yc = 0;  yd = 1;
BndryDescMat = [xa,yc,xb,yc, 0,-1; ... 
                xb,yc,xb,yd, 1, 0; ...
                xb,yd,xa,yd, 0, 1; ...
                xa,yd,xa,yc,-1, 0];
lambda = EqnBIC.lambda;
    mu = EqnBIC.mu;

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

%% Mesh generation and preparation 
% For a rectangular or smooth quadrilateral mesh 
n = 2^6;  nx = n;  ny= n;  delta = 0;  status = 2;
QuadriMesh = RectDom_QuadriMesh_GenLogiRectPtbSmth(...
  xa,xb,nx,yc,yd,ny,delta,status);
QuadriMesh = QuadriMesh_Enrich2(QuadriMesh,BndryDescMat);
QuadriMesh = QuadriMesh_Enrich4(QuadriMesh,BndryDescMat);

%% For time-marching: Now testing, just 1 step 
T = EqnBIC.T;  NT = n;  MT = NT;  Deltat = T/NT; 
tyme.T = T;  tyme.NT = NT;  tyme.MT = MT;  tyme.Deltat = Deltat;
tyme.tm = (0:NT)*Deltat;

%% Solving the coupled PDE system as a monolithic system (MS) 
[NumerDisp, NumerPres] ...
  = LinPoroElas_MS_QuadriIE_WGP02P02AC02P0_WGP0P0AC0_AsmSlv( ...
    EqnBIC,BndryCondTypeElas,BndryCondTypeDarcy,QuadriMesh,tyme,GAUSSQUAD);

%% Postprocessing to get quantities (at last step) for presentation 
% For Elas. 
[NumerDispEmAvgT, NumerStrsEmAvgT, NumerDivEmAvgT] ...
  = LinElas_WG_QuadriP02P02AC02P0_Presn( ... 
    EqnBIC.lambda, EqnBIC.mu, QuadriMesh, NumerDisp(:,MT+1), GAUSSQUAD);
% For Darcy 
NumerPresEmAvgT = NumerPres(1:QuadriMesh.NumEms, MT+1);

%% Computing errors when exact solutions are known 
[L2L2ErrDisp, L2L2ErrPres, L2L2ErrVel] ...
  = LinPoroElas_MS_QuadriIE_WGP02P02AC02P0_WGP0P0AC0_Err( ...
    EqnBIC,QuadriMesh,tyme,NumerDisp,NumerPres,GAUSSQUAD);

%% Presenting numerical and graphical results 
% [ min(NumerDivEmAvgT), max(NumerDivEmAvgT), ... 
%   min(NumerPresEmAvgT), max(NumerPresEmAvgT), ...
[L2L2ErrDisp, L2L2ErrPres, L2L2ErrVel]
show_QuadriMesh_ScaVecEm_mix(QuadriMesh, NumerDivEmAvgT, NumerDispEmAvgT, ...
  32, 'Numer. displacement & dilation at final time', 1);
show_QuadriMesh_ScaEm_ClrImg(QuadriMesh, NumerPresEmAvgT, ...
  34, 'Numer. pressure at final time');

return;