%% run_LinElas_CG_TriEP1.m 
% Solving lin.elasticity by 1st ord. enriched Lagrange elements on tri.mesh 
% 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;2;1;1];
% EqnBC = EqnBC_LinElas_YiM3AS2006_Ex2;
% 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 
n = 2^6;  nx = n;  ny = n;
TriMesh = RectDom_TriMesh_GenUnfmF(xa,xb,nx,yc,yd,ny,2);
TriMesh = TriMesh_Enrich2(TriMesh,BndryDescMat);
TriMesh = TriMesh_Enrich3(TriMesh,BndryDescMat);
TriMesh = TriMesh_Enrich4(TriMesh,BndryDescMat);
if TriMesh.NumEms<=2^6
  show_TriMesh_labels(TriMesh,11,0,0,0);
end

% %% Setting up the problem: Cook membrane 
% % Boundary condition type: 1 Dirichlet; 2 Neumann; 0 convenience 
% EqnBC = EqnBC_LinElas_CookMembrane;
% BndryCondType = [0;2;2;2;1];
% lambda = EqnBC.lambda;  
%     mu = EqnBC.mu;
% xa = 0;  xb = 48;  xc = 48;  xd = 0;
% ya = 0;  yb = 44;  yc = 60;  yd = 44;
% BndryDescMat = [ 0, 0, 48,44,   44/sqrt(48^2+44^2),-48/sqrt(48^2+44^2); ... 
%                 48,44, 48,60,                    1,                  0; ... 
%                 48,60,  0,44,  -16/sqrt(48^2+16^2), 48/sqrt(48^2+16^2); ... 
%                  0,44,  0, 0,                   -1,                 0];
% %% Mesh generation and preparation 
% n = 3 
% nx = 2^n;  ny = 2^(n-1);  delta = 0;  status = 2;
% TriMesh = UnitSqr_TriMesh_GenUnfm(n,n);
% TriMesh = TriMesh_Transform_TZdom(TriMesh,xa,ya,xb,yb,xc,yc,xd,yd);
% TriMesh = TriMesh_Enrich2(TriMesh,BndryDescMat);
% % TriMesh = TriMesh_Enrich4(TriMesh,BndryDescMat);
% if TriMesh.NumEms<=2^10 
%   show_TriMesh_labels(TriMesh,11,0,0,0);
% end

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

%% CG.EP1 for linear elasticity: Assembling & Solving
[sln] = LinElas_CG_TriEP1_AsmSlv(EqnBC,BndryCondType,TriMesh,GAUSSQUAD);

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

%% Computing errors when the exact solution is known 
[L2ErrDisp, L2ErrDiv, L2ErrStrs, L2RHS] ... 
  = LinElas_CG_TriEP1_Err(EqnBC,TriMesh,sln,GAUSSQUAD);
[AhErr] = LinElas_CG_TriEP1_ErrAh(EqnBC,TriMesh,sln,GAUSSQUAD);

%% Presenting numerical and graphical results
[AhErr, L2ErrDisp, L2ErrStrs, L2ErrDiv*lambda]  % /L2RHS
show_TriMesh_ScaVecEm_mix(TriMesh, NumerDivEmAvg, NumerDispEmAvg, ...
  23, 'CG.EP1 numer. displacement & divergence', 1);
show_TriMesh_ScaEm_ClrImg(TriMesh, NumerStrsEmAvg(:,1,2), ...
  24, 'CG.EP1 numer. shear stress \sigma_{xy}');

return;