%% run_Darcy_WG_QuadriQ0Q0RT0.m
% Solving Darcy eqn. by WG(Q0,Q0;RT[0])Quadri on a quadrilateral mesh 
% 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 (start/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_Ex12sinsin;
% BndryCondType = [0;2;2;2;1];
% EqnBC = EqnBC_Poisson_ArnoldBoffiFalkMCOM2002_Ex1;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_MuWangYeIJNAM2015_Case2;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_WheeYotovSINUM2006_Ex2;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_WihlerRivieJSC2011_Ex;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_ArboWheeYotovSINUM1997_Ex2;
% BndryCondType = [0;2;1;2;1];

%% Discretization: Mesh generation and preparation 
% 
% For a rectangular or smooth quadrilateral mesh 
n = 2^4;  nx = n;  ny= n;  delta = 1;  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 but asymp. paral. mesh (ArnoldBoffiFalk_MCOM_2002)
% m = 3  
% status = 3;
% QuadriMesh = RectDom_QuadriMesh_GenTZasymp(xa,xb,yc,yd,m,status);
% 
% % For a trapezoidal mesh in ArnoldBoffiFalk_MCOM_2002 
% n = 2^5;  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_Enrich3(QuadriMesh,BndryDescMat);
if QuadriMesh.NumEms<=2^8
  show_QuadriMesh_labels(QuadriMesh,11,0,0,0);
end

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

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

%% WG: QuadriQ0Q0RT0 quadri.mesh: Assembling and Solving 
% tic;
[sln,GlbMat,GlbRHS] = Darcy_WG_QuadriQ0Q0RT0_AsmSlv(...
  EqnBC,BndryCondType,QuadriMesh,PermK,GAUSSQUAD);
% toc

%% Computing numerical pressure, velocity, flux, ... 
[NumerPresEm, NumerVelCofRT0, NumerFlux, FluxDscp, LMCR] = ...
  Darcy_WG_QuadriQ0Q0RT0_PresVelFlux(QuadriMesh, PermK, sln);
NumerVelEmCntr = NumerVelCofRT0(:,1:2);

%% Computing errors when an exact solution is known 
[L2ErrPres,L2ErrVel,L2ErrFlux,L2ErrDiv] = Darcy_WG_QuadriQ0Q0RT0_Err( ...
  EqnBC,QuadriMesh,NumerPresEm,NumerVelCofRT0,GAUSSQUAD);

%% Presenting numerical and graphical results
[L2ErrPres, L2ErrVel, L2ErrFlux, L2ErrDiv]
% [L2ErrPres, L2ErrVel, L2ErrDiv]
% 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);
end

return;