%% run_NlinPoroElas_OS_QuadriIE_CGEQ1_WGP0P0AC0.m 
% Solving nonlinear poroelasticity via operator-splitting (OS) 
% 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) 
%    with dilation/stress-dependent permeability K(u) 
% First elas. then Darcy, since pressure changes faster 
% James Liu ColoState; 2012/07--2020/12 

% 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_sinsin;
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 = 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: Now testing, just 1 step 
T = 1e-3;  Deltat = 1e-3;  NT = T/Deltat;  MT = NT;
tyme.T = T;  tyme.NT = NT;  tyme.MT = MT;  tyme.Deltat = Deltat;
tyme.tm = (0:NT)*Deltat;

%% Solving the nonlinearly coupled PDE system via operator splitting (OS) 
% size(NumerDisp) = (2*NumNds+NumEgs) * MT 
% size(NumerPres) = (NumEms+NumEgs) * MT 
[NumerDisp, NumerPres] ... 
  = NlinPoroElas_OS_QuadriIE_CGBR1_WGP0P0AC0_AsmSlv( ... 
    EqnBIC, BndryCondTypeElas, BndryCondTypeDarcy, QuadriMesh, tyme, ... 
    GAUSSQUAD); 

%% Postprocessing: Extracting data 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 
show_QuadriMesh_ScaVecEm_mix(QuadriMesh, NumerDivEmAvgT, NumerDispEmCntrT, ...
  32, 'Numeri. displacement & dilation at final time', 1);
show_QuadriMesh_ScaEm_ClrImg(QuadriMesh, NumerPresEmAvgT, ...
  34, 'Numeri. pressure at final time');

return;