%% run_Darcy_WG_QuadriP1P0P02.m 
% Darcy: WG(P1,P0;P0^2)Quadri with stabilization 
% Ref.: Mu,Wang,Ye, JCAM, 2015 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2018/06 

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

%% Setting up the problem 
% A domain & its boundary: starting/ending points, outward unit normal vec.
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 type: 1 Dirichlet; 2 Neumann; 0 convenience 
EqnBC = EqnBC_Poisson_Ex11sinsin;
BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_WihlerRivieJSC2011_Ex;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_Ex50TestRT0;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_MuWangYeIJNAM2015_Case27
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_ArnoldBoffiFalkMCOM2002_Ex1;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_WheeYotovSINUM2006_Ex2;
% BndryCondType = [0;1;1;1;1];

%% Discretization: Mesh generation and preparation 
% % Rectangular or smooth quadrilateral meshes 
% n = 2^3;  nx = n;  ny= n;  delta = 0;  status = 2;
% QuadriMesh = RectDom_QuadriMesh_GenLogiRectPtbSmth(...
%   xa,xb,nx,yc,yd,ny,delta,status);
% 
% Trapezoidal meshes in ArnoldBoffiFalk_MCOM_2002 
n = 2^3;  nx = n;  ny = n;  alpha = 1/3;
QuadriMesh = RectDom_QuadriMesh_GenTrapezoid(xa,xb,nx,yc,yd,ny,alpha,2);
% 
% Asymptotically parallelogram trapezoidal mesh: ArnoldBoffiFalk_MCOM_2002
% m = 3;  status = 3;
% QuadriMesh = RectDom_QuadriMesh_GenTZasymp(xa,xb,yc,yd,m,status);
% 
QuadriMesh = QuadriMesh_Enrich2(QuadriMesh,BndryDescMat);
QuadriMesh = QuadriMesh_Enrich3(QuadriMesh,BndryDescMat);
if QuadriMesh.NumEms<=2^10
  show_QuadriMesh_labels(QuadriMesh,11,0,0,0);
end

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

%% Sampling the permeability tensor (2x2 matrix)
PermK = Darcy_SmplnPerm_QuadriMesh(EqnBC.fxnK,QuadriMesh,GAUSSQUAD);

% WG_QuadriP1P0P02t: Assembling and Solving 
% tic;
rho = 1;  % Stabilizer param. (DMP: discrete maximum principle) 
[sln,GlbMat,GlbRHS] = Darcy_WG_QuadriP1P0P02_AsmSlv( ...
  EqnBC,BndryCondType,QuadriMesh,PermK,GAUSSQUAD,rho);
% toc

%% Computing numerical pressure, velocity, flux, ... 
[NumerPresEm,NumerVelEmCntr,NumerFlux,LMCR,FluxDscp] = ...
  Darcy_WG_QuadriP1P0P02_PresVelFlux(QuadriMesh,PermK,sln);

%% Computing errors when the exact solution is known 
[L2ErrPres,L2ErrVel,ErrFlux] = ...
  Darcy_WG_QuadriP1P0P02_Err(EqnBC,QuadriMesh,PermK,sln,GAUSSQUAD);

%% Presenting numerical and graphical results
[L2ErrPres, min(NumerPresEm), max(NumerPresEm)]
% show_QuadriMesh_ScaEm_ClrImg(QuadriMesh, NumerPresEm, ...
%   21, 'WG Numerical pressure');
if QuadriMesh.NumEms<=2^14
  show_QuadriMesh_ScaVecEm_mix(QuadriMesh, NumerPresEm, NumerVelEmCntr,...
    23, 'WG Numerical pressure and velocity', 1);
  % show_scalar(FluxDscp, ...
  %   25, 'WG Numer. flux discrepancy across edges');
end

return;