%% run_LinElas_CG_TriEP1_Ldom.m (may also be called Gamma-domain) 
% Solving lin. elasticity on an L-shaped domain: (-1,1)^2\([0,1]x[-1,0]) 
% by CG.EP1 (1st ord. enriched Lagrange elements) on triangular meshes 
% James Liu, Graham Harper, ColoState; 2012/07--2020/11 

% 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_Alber3Comput2002_Ex1Gdom;
BndryCondType = [0;1;1;1;1;1;1];
lambda = EqnBC.lambda;  
    mu = EqnBC.mu;  
% The Gamma-shaped domain 
xa = -1;  xb = 1;  yc = -1;  yd = 1;
xm = 0.5*(xa+xb);  ym = 0.5*(yc+yd);
BndryDescMat = [xa,yc, xm,yc,  0, -1;
                xm,yc, xm,ym,  1,  0;
                xm,ym, xb,ym,  0, -1;
                xb,ym, xb,yd,  1,  0;
                xb,yd, xa,yd,  0,  1;
                xa,yd, xa,yc, -1,  0];

%% Mesh generation: First a quadrilateral mesh and then triangular meshes 
% Initial mesh: Mesh size h=1 
QuadriMesh = Ldom_QuadriMesh_GenInit();
QuadriMesh = QuadriMesh_Enrich1(QuadriMesh);
QuadriMesh = QuadriMesh_Enrich2(QuadriMesh,BndryDescMat);
if QuadriMesh.NumEms<=2^8
  show_QuadriMesh_labels(QuadriMesh,11,1,1,1);
end
% Refining the quadrilateral mesh m times so that mesh size h=1/(2^m) 
m = 2;
for k=1:m
  QuadriMesh = QuadriMesh_RefineReg(QuadriMesh);
  QuadriMesh = QuadriMesh_Enrich2(QuadriMesh,BndryDescMat);
end
QuadriMesh = QuadriMesh_Enrich3(QuadriMesh,BndryDescMat);
if QuadriMesh.NumEms<=2^8
  show_QuadriMesh_labels(QuadriMesh,12,1,1,1);
end
% Now generating a triangular mesh for the last quadrilateral mesh 
dir = 0;
TriMesh = TriMesh_RefinedFrom_QuadriMesh(QuadriMesh,dir);
% Enriching the triangular mesh 
TriMesh = TriMesh_Enrich2(TriMesh,BndryDescMat);
TriMesh = TriMesh_Enrich4(TriMesh,BndryDescMat);
if TriMesh.NumEms<=2^8
  show_TriMesh_labels(TriMesh,14,1,1,1);
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(:,2,2), ...
  24, 'CG.EP1 numer. shear stress \sigma_{yy}');

return;