%% run_NlinPoroElas_OS_QuadriIE_CGEQ1_WGP0P0AC0.m 
% Solving nonlinear poroelasticity via operator splitting (OS) 
% (Gauss-Seidel) on a quadrilateral mesh by implicit Euler with 
% -- CG Bernardi-Raugel(1) for elasticity (solid displacement)
% -- WG(P0,P0;AC0) for Darcy flow (fluid pressure) 
% A particular case: Permeability depending on dilation 
% James Liu ColoState; Zhuoran Wang, SYSU; 2012/07--2021/04 

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

%% Setting up the problem 
EqnBIC = EqnBIC_NlinPoroElas_LiuWangICCS2021_Ex3;
BndryCondTypeDarcy = [0;1;1;1;1];
BndryCondTypeElas  = [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(3,9,13);

%% Discretization: Mesh generation and preparation 
% For a rectangular or smooth quadrilateral mesh 
n = 2^5;  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);
% JL20171231: It seems QuadriMesh_Enrich3 is not needed for now  

%% 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 nlin. PDE system by operator-splitting (OS) / Gauss-Seidel (GS) 
tic
[NumerDisp,NumerPres,NumerVelCof] ... 
  = NlinPoroElas_OS_QuadriIE_CGEQ1_WGP0P0AC0_AsmSlv( ...
    EqnBIC,BndryCondTypeElas,BndryCondTypeDarcy,QuadriMesh,tyme,GAUSSQUAD);
toc;

%% Computing errors 
[L2ErrDspl_all_sqrt, L2ErrPres_all_sqrt, ... 
  H1ErrDspl_all, L2ErrVel_all_sqrt, L2ErrStrs_all_sqrt] ... 
  = LinPoroElas_OS_QuadriIE_CGEQ1_WGP0P0AC0_Err( ...
    EqnBIC,QuadriMesh,tyme,NumerDisp,NumerPres,NumerVelCof,GAUSSQUAD);

%% Computing physical quantities (at the last time step) for presentation 
[NumerDispEmT,NumerDivEmT,NumerStrsEmT] ...
  = LinElas_CG_QuadriEQ1_Presn(lambda,mu,QuadriMesh,NumerDisp(:,MT+1),GAUSSQUAD);
% Getting PermK in the last time step 
PermK = Darcy_SmplnNlinPerm_QuadriMesh(... 
  EqnBIC,QuadriMesh,GAUSSQUAD,NumerDisp(:,MT+1));
[NumerPresEmT,NumerVelEmT] ...
  = Darcy_WG_QuadriP0P0AC0_Presn( ... 
    QuadriMesh,PermK,NumerPres(:,MT+1),GAUSSQUAD);

%% Presenting numerical and graphical results 
% [min(NumerDivEmT), max(NumerDivEmT), min(NumerPresEmT), max(NumerPresEmT)];
[L2ErrDspl_all_sqrt, H1ErrDspl_all, L2ErrStrs_all_sqrt, ... 
  L2ErrPres_all_sqrt,L2ErrVel_all_sqrt]
show_QuadriMesh_ScaVecEm_mix(QuadriMesh, NumerPresEmT, NumerVelEmT, ...
  33, 'Numerical pressure and velocity at final time', 1);
% 
% figure(51);
% x = linspace(0,1,n);
% y = linspace(0,1,n);
% cv = [0.10 : 0.01 : 1.00];
% [C,h] = contour(x,y, reshape(NumerPresEmT,n,n), cv);
% clabel(C,h)
% xlim([0 1]);  
% ylim([0 1]);
% set(gca,'fontsize',18);
% title('Numer.pressure contours at final time','fontsize',18);
% axis square;

return;