%% run_Darcy_WG_HybridPQ0RT0.m
% Solving Darcy eqn. by WG(Q0,Q0;RT[0])Quadri+(P0,P0;RT0)Tri on hybrid 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 cond. type: 1 for Dirichlet; 2 for Neumann; 0 for convenience 
EqnBC = EqnBC_Poisson_Ex11sinsin;
BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_WheeYotovSINUM2006_Ex1;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_WheeYotovSINUM2006_Ex2;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_MuWangYeIJNAM2015_Case2;
% BndryCondType = [0;1;1;1;1];

%% Discretization: Mesh generation and preparation 
n = 2^3  % Choices: n=2^3,^4,^5,^6;  my=2,4,8,16 (a quarter of ny) 
nx = n;  ny = n;  my = (1/4)*ny;  delta = 0.5;  % Actually delta not used 
HybridMesh = RectDom_HybridMesh_GenLogiRectYlowPtb(xa,xb,nx,yc,yd,ny,my,delta);
HybridMesh = HybridMesh_Enrich2(HybridMesh,BndryDescMat);
if (HybridMesh.NumEms<=2^7)
  show_HybridMesh_labels(HybridMesh,11,0,0,0)
end

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

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

%% Assembling and Solving 
[sln,GlbMat,GlbRHS] = Darcy_WG_HybridPQ0RT0_AsmSlv(...
  EqnBC,BndryCondType,HybridMesh,PermK,GAUSSQUAD);

%% Computing numerical pressure, velocity, flux ... 
[NumerPresEm,NumerVelCofRT0,NumerFlux,LMCR,FluxDscp] = ...
  Darcy_WG_HybridPQ0RT0_PresVelFlux(HybridMesh,PermK,sln,GAUSSQUAD);
NumerVelEmCntr = [NumerVelCofRT0.quadri(:,1:2); NumerVelCofRT0.tri(:,1:2)];

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

%% Presenting numerical and graphical results 
[L2ErrPres, L2ErrVel, L2ErrFlux, min(NumerPresEm), max(NumerPresEm)]
% show_HybridMesh_ScaEm_ClrImg(HybridMesh, NumerPresEm,...
%   21, 'WG Numerical pressure');
show_HybridMesh_ScaVecEm_mix(HybridMesh, NumerPresEm, NumerVelEmCntr,... 
  23, 'WG Numerical pressure and velocity', 1);

return;