%% run_LinPoroElas_MS_TriIE_CGEP1_WGP0P0RT0_Mandel.m 
% Solving linear poroelasticity as a monolithic system (MS) 
% on a triangular mesh by implicit Euler with 
% -- CG Bernardi-Raugel(1) for linear elasticity (solid displacement) 
% -- WG(P0,P0;RT0) for Darcy flow (fluid pressure) 
% James Liu ColoState; Zhuoran Wang, SYSU; 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,
% 3 partial Dirichlet (symmetry in x-direction),
% 4 partial Dirichlet (symmetry in y-direction), i.e. u2 = 0,
% 5 partial Dirichlet (symmetry in y-direction): u2 = u2(y,t) 
EqnBIC = EqnBIC_LinPoroElas_Mandel;
BndryCondTypeElas  = [0,0; 0,4; 0,0; 2,5; 0,3];  % For bot./right/top/left 
BndryCondTypeDarcy = [0;2;1;2;2];
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;
     B = EqnBIC.B;
   nuu = EqnBIC.nuu;

%% Setting up Gaussian quadratures 
GAUSSQUAD = SetGaussQuad(3,9,13);

%% Mesh generation and preparation 
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_Enrich4(TriMesh,BndryDescMat);

%% For time-marching 
% T = EqnBIC.T;  NT = n;  Deltat = T/NT;  MT = NT;  
T = 50;  NT = 500;  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_TriIE_CGEP1_WGP0P0RT0_AsmSlv_Mandel( ...
    EqnBIC,BndryCondTypeElas,BndryCondTypeDarcy,TriMesh,tyme,GAUSSQUAD);

%% Computing ...
[NumerDispEmT,NumerStrsEmT,NumerDivEmT] ... 
  = LinElas_CG_TriBR1_Prsn(lambda,mu,TriMesh,NumerDisp(:,MT+1),GAUSSQUAD);
PermK = Darcy_SmplnPerm_TriMesh(EqnBIC.fxnK,TriMesh,GAUSSQUAD);
% t = 50
[NumerPresEmT1,NumerPresEmEgT1,NumerVelEmT1] = ...
  Darcy_WG_TriP0P0RT0_Prsn_mandel(TriMesh,PermK,NumerPres(:,MT+1),n);
NumerPresEmEgT1 = NumerPresEmEgT1/(B*(1+nuu)/3);
% t = 0.1
[NumerPresEmT2,NumerPresEmEgT2,NumerVelEmT2] = ...
  Darcy_WG_TriP0P0RT0_Prsn_mandel(TriMesh,PermK,NumerPres(:,2),n);
NumerPresEmEgT2 = NumerPresEmEgT2/(B*(1+nuu)/3);
% t = 1 
[NumerPresEmT3,NumerPresEmEgT3,NumerVelEmT3] = ...
  Darcy_WG_TriP0P0RT0_Prsn_mandel(TriMesh,PermK,NumerPres(:,11),n); 
NumerPresEmEgT3 = NumerPresEmEgT3/(B*(1+nuu)/3);
% t = 10
[NumerPresEmT4,NumerPresEmEgT4,NumerVelEmT4] = ...
  Darcy_WG_TriP0P0RT0_Prsn_mandel(TriMesh,PermK,NumerPres(:,101),n);
NumerPresEmEgT4 = NumerPresEmEgT4/(B*(1+nuu)/3);

%% Calculating nodal displacement values 
NumerDisp_ux_nodes = NumerDisp(1:2:2*TriMesh.NumNds,MT+1);
NumerDisp_uy_nodes = NumerDisp(2:2:2*TriMesh.NumNds,MT+1);

%% Presenting numerical and graphical results 
[min(NumerDivEmT),max(NumerDivEmT)]
[min(NumerPresEmEgT1),max(NumerPresEmEgT1)]
[min(NumerPresEmEgT2),max(NumerPresEmEgT2)]
% [min(NumerPresEmEgT3),max(NumerPresEmEgT3)]
% 
show_TriMesh_ScaVecEm_mix(TriMesh, NumerDivEmT, NumerDispEmT, ...
  23, 'Numerical displacement and dilation at final time', 1);
show_TriMesh_ScaVecEm_mix(TriMesh, NumerPresEmT1, NumerVelEmT1, ...
  33, 'Numerical pressure and velocity at final time', 1);

%% Plotting edge pressure (as a curve) 
% For t = 0.1
figure(1);
x = linspace(0,1,n+1);
num_pres_edge = NumerPresEmEgT2;
scatter(x, num_pres_edge, '*', 'black');
hold on;
% 
% For t = 1
figure(1);
x = linspace(0,1,n+1);
num_pres_edge = NumerPresEmEgT3;
scatter(x, num_pres_edge, '*', 'red');
% 
% For t = 10
figure(1);
x = linspace(0,1,n+1);
num_pres_edge = NumerPresEmEgT4;
scatter(x, num_pres_edge, '*', 'blue');
% 
% For t = 50
figure(1);
x = linspace(0,1,n+1);
num_pres_edge = NumerPresEmEgT1;
scatter(x, num_pres_edge, '*', 'green');
hold on;

fig = legend(... 
  'Anal.sln. at T=0.1', ... 
  'Anal.sln. at T=1', ... 
  'Anal.sln. at T=10', ...
  'Anal.sln. at T=50', ...
  'Numer.sln at T=0.1', ... 
  'Numer.sln at T=1', ... 
  'Numer.sln at T=10', ... 
  'Numer.sln at T=50', ... 
  'Location', 'southwest');
set(fig, 'FontSize', 6);  
xlabel('x');  ylabel('Numer.pressure');

return;