%% run_LinElas_WG_QuadriP02P02AC02P0.m 
% Solving lin. elasticity by WG(P0^2,P0^2;AC0^2,P0) on a quadri. mesh
%   Discrete weak gradient:   (P0^2,P0^2) -> AC0^2 
%   Discrete weak divergence: (P0^2,P0^2) -> P0 
% James Liu, ColoState; 2017/02--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_CarsScheIMAJNA2015_Ex1;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_LinElas_HansboLarsonM2AN2003_ElasCavity;
% 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 normal) 
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;  status = 2;
QuadriMesh = RectDom_QuadriMesh_GenLogiRectPtbSmth(...
  xa,xb,nx,yc,yd,ny,delta,status);
% % For a randomly perturbed rough quadrilateral mesh 
% n = 2^3;  nx = n;  ny = n;  delta = 0.25;  status = 2;
% QuadriMesh = RectDom_QuadriMesh_GenLogiRectPtbRand(...
%   xa,xb,nx,yc,yd,ny,delta,status);
% % For a trapezoidal mesh in ArnoldBoffiFalk_MCOM_2002 
% n = 2^7;  nx = n;  ny= n;  alpha = 0.25;
% QuadriMesh = RectDom_QuadriMesh_GenTrapezoid(xa,xb,nx,yc,yd,ny,alpha,2);
QuadriMesh = QuadriMesh_Enrich2(QuadriMesh,BndryDescMat);
QuadriMesh = QuadriMesh_Enrich3(QuadriMesh,BndryDescMat);
% JL20160719: It seems QuadriMesh_Enrich3 is not needed for now  
if QuadriMesh.NumEms<=2^6; show_QuadriMesh_labels(QuadriMesh,11,0,0,0); end

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

%% WG(P0^2,P0^2;AC02,P0)Quadri for linear elasticity: Assembling & Solving
[sln] = LinElas_WG_QuadriP02P02AC02P0_AsmSlv( ...
  EqnBC,BndryCondType,QuadriMesh,GAUSSQUAD);

% %% Computing numer. displacement, stress, divergence: Coeffs. in resp.bas. 
% [NumerDispCofP02,NumerStrnCof,NumerDivCofP0] = ...
%   LinElas_WG_QuadriP02P02AC02P0_DispDivStrs(QuadriMesh,sln,GAUSSQUAD);
% NumerDispEmAvg = NumerDispCofP02;
% NumerDivEmAvg = NumerDivCofP0;
% NumerStrsEmAvg = zeros(QuadriMesh.NumEms,2,2);
% Wavgc = Hdiv_QuadriAC02_NmlzPiolaBas_FxnValAvg(QuadriMesh,0.5,0.5);
% NumerStrsEmAvg(:,1,1) = ... 
%   (2*mu)*sum(NumerStrnCof.*squeeze(Wavgc(:,:,1,1)),2) + lambda*NumerDivCofP0;
% NumerStrsEmAvg(:,2,2) = ... 
%   (2*mu)*sum(NumerStrnCof.*squeeze(Wavgc(:,:,2,2)),2) + lambda*NumerDivCofP0;
% NumerStrsEmAvg(:,1,2) = (2*mu)*sum(NumerStrnCof.*squeeze(Wavgc(:,:,1,2)),2);
% NumerStrsEmAvg(:,2,1) =  NumerStrsEmAvg(:,1,2);

%% Computing quantities for presenttaion 
[NumerDispEmAvg,NumerStrsEmAvg,NumerDivEmAvg] ...
  = LinElas_WG_QuadriP02P02AC02P0_Presn(lambda,mu,QuadriMesh,sln,GAUSSQUAD);

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

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

return;