%% run_LinPoroElas_MS_TriIE_CGEP1_WGP0P0RT0.m 
% Solving linear poroelasticity as a monolithic system (MS) 
% on a triangular mesh with implicit Euler (IE) and 
% -- CG.EP1 (enriched Lagrange P1) for lin.elasticity (solid displacement)
% -- WG(P0,P0;RT0) for Darcy flow (fluid pressure) 
% 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_LiuWangICCS2021_Ex0sinsinsin;
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;  
TriMesh = RectDom_TriMesh_GenUnfm(xa,xb,nx,yc,yd,ny,2);
TriMesh = TriMesh_Enrich2(TriMesh,BndryDescMat);
TriMesh = TriMesh_Enrich3(TriMesh,BndryDescMat);
TriMesh = TriMesh_Enrich4(TriMesh,BndryDescMat);
% JL20201226: It seems TriMesh_Enrich3 is not needed for now 
if TriMesh.NumEms<=2^6
  show_TriMesh_labels(TriMesh,11,0,0,0);
end

%% For time-marching: 
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) 
% size(NumerDisp) = (2*NumNds+NumEgs) * NT 
% size(NumerPres) = (NumEms+NumEgs) * NT 
[NumerDisp,NumerPres]... 
  = LinPoroElas_MS_TriIE_CGEP1_WGP0P0RT0_AsmSlv(... 
    EqnBIC,BndryCondTypeElas,BndryCondTypeDarcy,TriMesh,tyme,GAUSSQUAD);
%% Computing errors when exact solutions are known 
[L2L2ErrDisp, L2L2ErrPres] ... 
  = LinPoroElas_MS_TriIE_CGEP1_WGP0P0RT0_Err( ... 
    EqnBIC,TriMesh,tyme,NumerDisp,NumerPres,GAUSSQUAD);

%% Postprocessing: Extracting data for presentation 
[NumerDispEmCntrT, NumerDivEmAvgT, NumerStrsEmAvgT] ... 
  = LinElas_CG_TriEP1_Presn(lambda,mu,TriMesh, ... 
    NumerDisp(:,MT+1),GAUSSQUAD); 
NumerPresEmAvgT = NumerPres(1:TriMesh.NumEms, MT+1);

%% Presenting numerical and graphical results 
[min(NumerDivEmAvgT), max(NumerDivEmAvgT), ... 
  min(NumerPresEmAvgT), max(NumerPresEmAvgT), ... 
  L2L2ErrDisp, L2L2ErrPres] 
show_TriMesh_ScaVecEm_mix(TriMesh, NumerDivEmAvgT, NumerDispEmCntrT, ...
  32, 'CG.EP1 Numer.displacement & dilation at final time', 1);
show_TriMesh_ScaEm_ClrImg(TriMesh, NumerPresEmAvgT, ...
  34, 'WG Numer.pressure at final time');

return;