%% run_LinPoroElas_MS_QuadriIE_CGEQ1_WGP0P0AC0.m 
% Solving linear poroelasticity as a monolithic system (MS) 
% on a quadrilateral mesh with implicit Euler (IE) and 
% -- CG.EQ1 (enriched Lagrange Q1) for lin.elasticity (solid displacement)
% -- WG(P0,P0;AC0) for Darcy flow (fluid pressure) 
% First elas. then Darcy, since pressure changes faster 
% James Liu ColoState; 2017/07--2021/04 

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

%% Setting up the problem 
% Boundary condition types: 0 convenience, 1 Dirichlet, 2 Neumann 
% EqnBIC = EqnBIC_LinPoroElas_KansRivieJSC2018_ExVar;
% BndryCondTypeElas  = [0;1;1;1;1];
% BndryCondTypeDarcy = [0;1;1;1;1];
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 preparation 
% For a rectangular or smooth quadrilateral mesh 
n = 2^4;  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_CGEQ1_WGP0P0AC0_AsmSlv(... 
    EqnBIC,BndryCondTypeElas,BndryCondTypeDarcy,QuadriMesh,tyme,GAUSSQUAD); 

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

%% Postprocessing: Obtaining quantities (at last step) for presentation 
[NumerDispEmCntrT, NumerStrsEmAvgT, NumerDivEmAvgT] ... 
  = LinElas_CG_QuadriEQ1_Presn(lambda, mu, QuadriMesh, ... 
    NumerDisp(:,MT+1), GAUSSQUAD); 
NumerPresEmAvgT = NumerPres(1:QuadriMesh.NumEms, MT+1);

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

return;