%% run_LinPoroElas_MS_TriIE_CGEP1_WGP0P0RT0_mandel.m 
% Solving linear poroelasticity as a monolithic system (MS)
% on a triangular mesh by implicit Euler 
% -- CG Bernardi-Raugel(1) for elasticity eqn. (solid displacement)
% -- WG(P0,P0;RT0) for Darcy eqn. (fluid pressure) 
% James Liu ColoState; 2012/07--2020/07 

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

%% Setting up the problem 
% Setting up the domain 
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];
% Boundary condition types: 0 convenience, 1 Dirichlet, 2 Neumann,
% 3 partial Dirichlet (symmetry in x-direction),
% 4 partial Dirichlet (symmetry in y-direction), and u2 = 0,
% 5 partial Dirichlet (symmetry in y-direction), and u2 = u2(y,t),

EqnBIC = EqnBIC_LinPoroElas_Mandel;
BndryCondTypeDarcy = [0;2;1;2;2];
% BndryCondTypeElas  = [0;4;2;5;3];  
BndryCondTypeElas  = [0,0;0,4;0,0;2,5;0,3];  
% BndryCondTypeElas  = [0,0;0,4;0,0;0,5;0,3]; 

% EqnBIC = EqnBIC_LinPoroElas_sinsin;
% BndryCondTypeDarcy = [0;1;1;1;1];
% BndryCondTypeElas  = [0,0;1,0;1,0;1,0;1,0];
%
lambda = EqnBIC.lambda;
    mu = EqnBIC.mu;
     B = EqnBIC.B;
   nuu = EqnBIC.nuu;


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

%% Discretization: Mesh generation and preparation 
% For a rectangular or smooth quadrilateral mesh 
% m = 2;
% n = 2^5;  
% nx = n;  ny= n;  delta = 0;  status = 2;
% TriMesh = RectDom_TriMesh_GenLogiRectPtbSmth(...
%   xa,xb,nx,yc,yd,ny,delta,status);
% % 
% % For a trapezoidal but asymp. paral. mesh (ArnoldBoffiFalk_MCOM_2002)
% % m = 6;  n = 2^m;  status = 3;
% % TriMesh = RectDom_TriMesh_GenTZasymp(xa,xb,yc,yd,m,status);
% % 
% % For a trapezoidal mesh in ArnoldBoffiFalk_MCOM_2002 
% % n = 2^5;  nx = n;  ny= n;  alpha = 0;
% % TriMesh = RectDom_TriMesh_GenTrapezoid(xa,xb,nx,yc,yd,ny,alpha,2);
% 
% TriMesh = TriMesh_Enrich2(TriMesh,BndryDescMat);
% TriMesh = TriMesh_Enrich4(TriMesh,BndryDescMat);
% % JL20171231: It seems TriMesh_Enrich3 is not needed for now  
% if TriMesh.NumEms<=2^7
%   show_TriMesh_labels(TriMesh,11,1,1,1);
% end

m = 3;
n = 2^m;  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);
if TriMesh.NumEms<=2^7
  show_TriMesh_labels(TriMesh,11,1,1,1);
end

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


%% Computing errors 
% [L2ErrDspl, L2ErrPres, H1ErrDspl] = LinPoroElas_MS_QuadriIE_CGBR1_WGP0P0AC0_Err( ...
%   EqnBIC,TriMesh,tyme,NumerDspl,NumerPres,GAUSSQUAD);


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

%% To plot nodal displacement values.
NumerDspl_ux_nodes = NumerDspl(1:2:2*TriMesh.NumNds,MT+1);
NumerDspl_uy_nodes = NumerDspl(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)]
% 
% [L2ErrDspl, L2ErrPres, H1ErrDspl]
show_TriMesh_ScaVecEm_mix(TriMesh, NumerDivEmT, NumerDsplEmT, ...
  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);
%
% show_TriMesh_ScaEm_ClrImg(TriMesh,NumerStrsEmT(:,2,2), ...
%   43, 'Numerical effective stress(yy) at final time');
% show_TriMesh_ScaEm_ClrImg(TriMesh,NumerStrsEmT(:,1,1), ...
%   13, 'Numerical effective stress(xx) at final time');
%
% print -f23 -dpng ExSandwichedLowPerm_RectIE_BR1WG0_DsplDiv_n32t.png 
% print -f33 -dpng ExSandwichedLowPerm_RectIE_BR1WG0_PresVel_n32t.png 
% % hold on 
% figure(51);
% x = linspace(0,1,n);
% y = linspace(0,1,n);
% % cv = [-3.0 : 0.1 : 3.0];
% % cv = [[0.00 : 0.20: 0.90], [0.90 : 0.01 : 1.00]];
% cv = [0.10 : 0.01: 1.00];
% [C,h] = contour(x,y,reshape(NumerPresEmT1,n,n),cv);  % /(B*F*(1+nuu)/3)???
% clabel(C,h)
% xlim([0 1]);  ylim([0 1]);
% set(gca,'fontsize',18);
% title('Numerical pressure contours at T=','fontsize',18);
% axis square;
% colormap cool;
% print -f51 -dpng ExSandwichedLowPerm_RectIE_BR1WG0_PresCntur_n32t.png


% %% displacement contour
% figure(52);
% x = linspace(0,1,n+1);
% y = linspace(0,1,n+1);
% % cv = [-3.0 : 0.1 : 3.0];
% % cv = [[0.00 : 0.20: 0.90], [0.90 : 0.01 : 1.00]];
% [C,h] = contour(x,y,reshape(NumerDspl_ux_nodes,n+1,n+1),10);  % /(B*F*(1+nuu)/3)???
% clabel(C,h)
% xlim([0 1]);  ylim([0 1]);
% set(gca,'fontsize',18);
% title('Numerical ux contours at T=','fontsize',18);
% axis square;
% 
% figure(53);
% x = linspace(0,1,n+1);
% y = linspace(0,1,n+1);
% % cv = [-3.0 : 0.1 : 3.0];
% % cv = [[0.00 : 0.20: 0.90], [0.90 : 0.01 : 1.00]];
% [C,h] = contour(x,y,reshape(NumerDspl_uy_nodes,n+1,n+1),10);  % /(B*F*(1+nuu)/3)???
% clabel(C,h)
% xlim([0 1]);  ylim([0 1]);
% set(gca,'fontsize',18);
% title('Numerical uy contours at T=','fontsize',18);
% axis square;

% %% velocity contour
% figure(54);
% x = linspace(0,1,n);
% y = linspace(0,1,n);
% % cv = [-3.0 : 0.1 : 3.0];
% % cv = [[0.00 : 0.20: 0.90], [0.90 : 0.01 : 1.00]];
% [C,h] = contour(x,y,reshape(NumerVelEmT1(:,1),n,n),10);  % /(B*F*(1+nuu)/3)???
% clabel(C,h)
% xlim([0 1]);  ylim([0 1]);
% set(gca,'fontsize',18);
% title('Numerical velocity_x contours at T=','fontsize',18);
% axis square;
% 
% 
% figure(55);
% x = linspace(0,1,n);
% y = linspace(0,1,n);
% % cv = [-3.0 : 0.1 : 3.0];
% % cv = [[0.00 : 0.20: 0.90], [0.90 : 0.01 : 1.00]];
% [C,h] = contour(x,y,reshape(NumerVelEmT1(:,2),n,n),10);  % /(B*F*(1+nuu)/3)???
% clabel(C,h)
% xlim([0 1]);  ylim([0 1]);
% set(gca,'fontsize',18);
% title('Numerical velocity_y contours at T=','fontsize',18);
% axis square;
% 
% %% sigma contour
% figure(56);
% x = linspace(0,1,n);
% y = linspace(0,1,n);
% % cv = [-3.0 : 0.1 : 3.0];
% % cv = [[0.00 : 0.20: 0.90], [0.90 : 0.01 : 1.00]];
% [C,h] = contour(x,y,reshape(NumerStrsEmT(:,2,2),n,n),10);  % /(B*F*(1+nuu)/3)???
% clabel(C,h)
% xlim([0 1]);  ylim([0 1]);
% set(gca,'fontsize',18);
% title('Numerical sigma_yy contours at T=','fontsize',18);
% axis square;
% 
% figure(57);
% x = linspace(0,1,n);
% y = linspace(0,1,n);
% % cv = [-3.0 : 0.1 : 3.0];
% % cv = [[0.00 : 0.20: 0.90], [0.90 : 0.01 : 1.00]];
% [C,h] = contour(x,y,reshape(NumerStrsEmT(:,1,1),n,n),10);  % /(B*F*(1+nuu)/3)???
% clabel(C,h)
% xlim([0 1]);  ylim([0 1]);
% set(gca,'fontsize',18);
% title('Numerical sigma_xx contours at T=','fontsize',18);
% axis square;


%% pressure edge contour
% t = 0.1
figure(1);
x = linspace(0,1,n+1);
num_pres_edge = NumerPresEmEgT2;
% num_pres_center = num_pres_center;
scatter(x,num_pres_edge,'*','black');
hold on;

% t = 1
figure(1);
x = linspace(0,1,n+1);
num_pres_edge = NumerPresEmEgT3;
% num_pres_center = num_pres_center;
scatter(x,num_pres_edge,'*','red');

% t = 10
figure(1);
x = linspace(0,1,n+1);
num_pres_edge = NumerPresEmEgT4;
% num_pres_center = num_pres_center;
scatter(x,num_pres_edge,'*','blue');

% 
% t = 50
figure(1);
% x = linspace(1/n,(n-1)/n,n);
% num_pres_center = NumerPresEmT1(1:2^m:4^m);
x = linspace(0,1,n+1);
num_pres_edge = NumerPresEmEgT1;
% num_pres_edge = num_pres_edge;
scatter(x,num_pres_edge,'*','green');
hold on;


fig = legend('Analytical soln. at T = 0.1','Analytical soln. at T = 1','Analytical soln. at T = 10','Analytical soln. at T = 50', ...
       'Num. soln at T = 0.1','Num. soln at T = 1','Num. soln at T = 10','Num. soln at T = 50','Location','southwest');
% fig = legend('Analytical soln. at T = 0.01','Analytical soln. at T = 0.1','Analytical soln. at T = 1','Analytical soln. at T = 5', ...
%        'Num. soln at T = 0.01','Num. soln at T = 0.1','Num. soln at T = 1','Num. soln at T = 5','Location','northeast');
set(fig,'FontSize',6);  
xlabel('x');
ylabel('Num. pressure');


return;