%% run_LinElas_CG_QuadriEQ1.m 
% Solving lin. elasticity on a quadri. mesh 
% by enriched Lagrangian Q1^2 elements (with "reduced integration") 
% James Liu, ColoState; 2012/07--2021/02 

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

%% Setting up the problem to be solved 
% Boundary condition types: 1 Dirichlet; 2 Neumann, 0 convenience 
EqnBC = EqnBC_LinElas_BrenSINUM1993_Ex;
BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_LinElas_Liu_Ex2;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_LinElas_Wang3ZhangJCAM2016_Test3Locking;
% BndryCondType = [0;1;1;1;1];
lambda = EqnBC.lambda; 
    mu = EqnBC.mu;
% The domain & its boundary (starting/ending points, outward unit normals) 
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];

%% Mesh generation and preparation 
% % For a rectangular or smooth quadrilateral mesh 
n = 2^6;  nx = n;  ny= n;  delta = 0.0;  status = 2;
QuadriMesh = RectDom_QuadriMesh_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;
% QuadriMesh = RectDom_QuadriMesh_GenTZasymp(xa,xb,yc,yd,m,status);
% 
% % For a trapezoidal mesh in ArnoldBoffiFalk_MCOM_2002 
% n = 2^6;  nx = n;  ny= n;  alpha = 1/3;
% QuadriMesh = RectDom_QuadriMesh_GenTrapezoid(xa,xb,nx,yc,yd,ny,alpha,2);
% 
QuadriMesh = QuadriMesh_Enrich2(QuadriMesh,BndryDescMat);
QuadriMesh = QuadriMesh_Enrich4(QuadriMesh,BndryDescMat);
% JL20160719: It seems QuadriMesh_Enrich3 is not needed for now  
if QuadriMesh.NumEms<=2^6 
  show_QuadriMesh_labels(QuadriMesh,11,0,1,0);
end

%% Setting up Gaussian quadratures 
GAUSSQUAD = SetGaussQuad(5,25,13);

%% CG.EQ1 w/ red.integration for lin.elasticity: Assembling & Solving
[sln] = LinElas_CG_QuadriEQ1_AsmSlv( ...
  EqnBC,BndryCondType,QuadriMesh,GAUSSQUAD);

%% Comput. numeri. displacement, stress, dilation em.avg. for presentation 
[NumerDispEmAvg,NumerStrsEmAvg,NumerDivEmAvg] = ...
  LinElas_CG_QuadriEQ1_Presn(lambda,mu,QuadriMesh,sln,GAUSSQUAD);

%% Computing errors when an exact solution is known 
[L2ErrDisp,L2ErrStrs,L2ErrDiv,L2RHS] = ...
  LinElas_CG_QuadriEQ1_Err(EqnBC,QuadriMesh,sln,GAUSSQUAD);

%% Presenting numerical and graphical results 
[L2ErrDisp,L2ErrStrs,lambda*L2ErrDiv]/L2RHS 
show_QuadriMesh_ScaVecEm_mix(QuadriMesh, NumerDivEmAvg, NumerDispEmAvg, ...
  23, 'CG.EQ1 numerical displacement & dilation', 1);
show_QuadriMesh_ScaEm_ClrImg(QuadriMesh,NumerStrsEmAvg(:,1,1), ...
  24, 'CG.EQ1 numerical stress \sigma_{xx}');

return;